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ABSTRACT 


This  investigation,  develops  an  axisymmetric  heat  transfer 
combustion  model  of  a  porous  medium  within  a  circular  cylin¬ 
der.  System  flow  is  governed  by  Darcy's  law.  Carbon  and 
air  properties  are  treated  as  variables  of  temperature.  A 
combined  cont inuity-Darcy  equation,  an  oxygen  mass  balance 
equation,  and  energy  balance  equations  (one  each)  for  air 
and  carbon,  describe  the  conservation  laws  of  the  system. 
Transport  mechanisms  for  oxygen  mass  transfer  are  molecular 
diffusion  and  convective  transport,  and  an  oxygen  consumption 
term  to  account  for  combustion  is  included.  Heat  transfer 
mechanisms  included  in  the  model  are  conduction  and  convec¬ 
tion.  Radiation  is  accounted  for  at  applicable  boundaries 
only.  Nonvolatile  combustion  is  accounted  for  in  the  carbon 
energy  and  oxygen  mass  balance  equations  as  a  heat  generation 
term  of  Arrhenius  type.  The  numerical  solution  of  four 
coupled,  nonlinear,  transient  partial  differential  field 
equations  is  accomplished  using  the  Galerkin  formulation  of 
the  Finite  Element  Method.  The  effect  of  porosity  on  system 
behavior  is  examined.  ^  ' 
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increase  with  increasing  Reynolds  numbers.  Boffa  [Ref.  12] 
has  shown  that  for  a  fixed  Reynolds  number,  inertial  effects 


diminish  with  increasing  air  temperature.  Darcy's  law  for 
two-dimens ional  flow  is, 


■m 


(dP 

ldr 


r  +( 


dP 

dz 


"  P 


g)  z , 


(3.5) 


where  Q  is  the  filter  velocity  or  the  volumetric  flow  rate 
per  unit  cross-sectional  area;  m  is  the  specific  permeability 

w.  -\ 

g  is  the  gravitational  acceleration;  r  and  z  are  the 
unit  vectors  in  the  r  and  z  directions,  respectively;  and 
dP/dr  and  dP/dz  are  the  pressure  gradients  in  the  r  and  z 
directions,  respectively.  The  assumption  here  is  that  the  r 
(radial)  and  z  (axial)  velocity  components  each  react  to  the 
pressure  independent  of  the  other.  The  specific  permeability 
of  the  porous  medium  used  in  the  model  is, 


m 


96  V 


(3.6) 


Expression  3.6  is  based  on  a  cappilaric-serial  model  given 
by  Scheidegger  (Ref.  10].  Expressions  for  permeability  vary 
with  the  physical  assumptions  of  the  flow  paths  incorporated 
into  the  model.  Permeability  also  varies  with  tne  oore  size 
distribution  assumed.  Physically,  permeability  and  porosity 
are  not  related  [Ref.  9].  Porosity  is  a  quantifiable 
property  of  the  porous  medium,  whereas  permeability  is  a 
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as  0.99.  Changes  in  particle  diameter  are  incorporated  in 
the  model  and  are  discussed  in  Section  III.D  with  carbon 
combustion . 

B.  DARCY’S  LAW  AND  PORE  VELOCITY 

The  Reynolds  number  for  porous  media  is  defined  by 

a  s  d 

Re  =  — -  (3.4) 

U 

where  s  is  the  local  pore  velocity,  p  is  the  mass  density  of 
air,  and  „  is  the  dynamic  viscosity.  The  magnitude  of  the 
Reynolds  number  indicates  whether  fluid  motion  is  dominated 
by  molecular,  viscous,  or  inertial  effects.  Most  investi¬ 
gations  of  flow  through  porous  media  indicate  flow  regimes 
of  viscous  and  inertially  dominated  flows.  The  Navier- 
Stok.es  equations  apply  to  fluid  motion  possessing  such 
Reynolds  numbers.  Because  of  the  geometry  involved  in  a 
consolidated  (rigid)  porous  medium  and  the  no-slip  boundary 
condition  (i.e.,  s  =  0  at  a  solid-fluid  interface),  solution 
of  the  Navier-Stokes  equations  is  difficult  for  a  porous 
medium.  Scheidegger  [Ref.  10]  points  out  that  extensive 
experimental  work  with  porous  media  indicates  fluid  flow  is 
governed  by  Darcy's  law  for  the  range  of  Reynolds  numbers 
where  viscous  effects  dominate.  The  upper  limit  (velocity) 
Reynolds  number  in  these  investigations  is  subject  to  dis¬ 
agreement  and  varies  from  0.1  to  75. 
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Inertial  effects 


for  spherical  particles.  Equation  3.3  is  sometimes  used 
and  assumes  one-half  of  the  total  internal  surface  area 


\ 


I 

I 


is  effective  for  convective  heat  transfer.  The  fractional 
amount  of  total  area  is  an  estimate  based  on  Fontenot's  [Ref.  37] 
experimental  results  and  does  not  generally  apply  to  porous 
media  [Ref.  9],  The  Kozeny  relations,  alternate  expressions 
of  the  specific  internal  area,  are  discussed  by  Scheidegger 
[Ref.  10].  Advantages  of  the  Kozeny  relations  are  fair 
agreement  with  experimental  values  and  calculations  that  are 
independent  of  particle  shape.  A  disadvantage  is  the 
failure  of  the  relations  to  predict  accurate  values  of  z  at 
high  values  of  porosity  [Ref.  9].  For  the  geometric 
configuration  of  Figure  3.1,  the  tortuosity  depends  on  the 
ratio  d/D.  Carman  [Ref.  11]  presents  a  table  of  measured 
tortuosity  factors  of  various  materials  and  geometries,  and 
points  out  differences  between  analytical  determinations  of 
tortuosity.  In  this  study,  his  recommended  value  of  1.4  is 
used.  Particle  size  decreases  with  consumption.  As  a 
result,  thermophysical  properties  which  depend  on  particle 
diameter,  as  well  as  temperature,  are  functions  of  time  and 
space.  The  numerical  model  presented  assumes  matrix 
rigidity  as  particle  diameter  decreases.  This  assumption 
gives  rise  to  an  increase  in  porosity  during  combustion. 

Although  Scheidegger  [Ref.  10]  in  his  discussion  of  the 
packing  theory  of  spheres,  reports  .875  porosity  as  the 
threshold  for  stability  in  a  porous  matrix,  Carman  [Ref.  11] 
reports  on  investigations  performed  on  porosities  as  high 
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path  to  the  straight  (line)  path  displacement  of  the 
particle.  Permeability,  a  measure  of  hydraulic  conductivity, 
is  a  property  of  the  porous  medium  that  depends  upon  the 
four  characteristics  mentioned  above.  Scheidegger  [Ref.  10] 
presents  methods  used  to  measure  the  properties  of  porous 
media.  Methods  discussed  are  essentially  experimental  in 
nature . 

The  porous  medium  was  modelled  as  shown  in  Figure  3.1. 

In  Figure  3.1,  D  is  the  particle  center- to-center  distance, 
and  d  is  the  particle  diameter.  From  the  idealized  geometry, 
the  porosity  for  spherical  particles  is, 

P  “  1  -  f(i)3  (3-15 

The  pore  diameter  is  obtained  from  an  expression  proposed 
by  Carman  [Ref.  11], 


Equation  3.2  is  analogous  to  a  more  familiar  form  of  mean 
hydraulic  diameter,  4v/A,  where  v  is  the  void  volume  and  A 
is  the  wetted  surface  area.  The  specific  internal  area,  Z, 
based  on  the  idealized  geometry  of  Figure  3.1  may  be  expressed 
as , 


Z 


(3.3) 
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III. 


THEORY  AND  BACKGROUND 


A.  DESCRIPTION  OF  THE  POROUS  MEDIUM 

In  this  work,  a  porous  medium  is  considered  to  be  a 
solid  containing  interconnecting  pores  that  allow  fluid  to 
permeate  and  flow  through  the  solid.  Either  of  two  classes 
of  porous  media,  consolidated  (solid  and  rigid)  and  uncon¬ 
solidated  (comprised  of  discrete  particles  as  found  in 
granular  beds)  may  be  considered.  Each  class  of  porous 
media  may  have  isotropic  nonhomogeneous  properties. 

The  common  characteristics  of  all  porous  media  are: 

(1)  porosity,  (2)  specific  internal  area,  (3)  pore  diameter, 
and  (4)  tortuosity.  Porosity,  p,  is  defined  as  the  ratio  of 
void  volume  to  total  volume.  n’he  specific  internal  area,  z, 
is  the  ratio  of  internal  surface  area  to  bulk  volume.  In 
general,  the  distribution  of  pore  size  in  a  porous  medium  is 
random  (nonhomogeneous)  and  dynamic  (subject  to  small 
strains  induced  by  the  transient  pressure  field) .  The  situa¬ 
tion  motivates  the  investigator  to  treat  the  porous  medium 
as  a  continuum  possessing  idealistic  geometrical  properties 
of  porosity  and  pore  diameter.  The  specific  internal  area 
is  obtained  from  the  porosity  model.  Though  many  conventions 
exist  for  the  definition  of  a  conceptual  pore  diameter,  in 
this  study  a  hydraulic  diameter  analogy  is  employed.  The 
tortuosity,  t,  of  a  porous  medium  is  the  ratio  of  the  flow 
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transfer  mechanism  of:  (1)  molecular  diffusion,  and 

(2)  convective  transport  of  species.  An  oxygen  consumption 

term  due  to  combustion  is  included. 

The  fourth  field  equation  involves  the  system  pressure 
gradient.  The  equation  is  a  combined  Darcy's  law  and  air 
mass  continuity  equation  for  the  system. 

The  conservation  equations  describing  the  system  field 
are  four  transient,  coupled,  nonlinear  partial  differential 
equations.  The  four  equations  are  solved  by  a  two-dimensional 
Galerkin  formulation  of  the  Finite  Element  Method.  The  inte¬ 
gration  scheme  used  for  this  highly  stiff  system  is  one 
presented  by  Gear  and  modified  by  Franke  [Ref.  30].  The 
scheme  is  especially  suited  to  systems  of  implicit  and  stiff 
differential  equations.  The  integration  scheme  is  used  in 
conjunction  with  an  optimal  compact  storage  scheme  proposed 
by  Franke  and  Salinas  [Ref.  29]  . 


PROBLEM  DESCRIPTION 


II  . 

The  problem  under  investigation  is  that  of  combustion 
of  a  porous  fuel  (carbon)  imbedded  in  a  cylindrical  container. 
Pressure  gradients  induce  convective  currents  through  the 
medium.  The  air  flow  produces  two  opposing  effects: 

(1)  internal  convective  heat  transfer,  and  (2)  a  supply  of 
oxygen  to  promote  heat  generation  through  combustion. 
Extinguishment  or  sustained  combustion  of  the  porous  medium 
depend  on  the  interaction  of  these  effects.  If  heat  transfer 
dominates,  the  combustion  will  proceed  to  extinguishment, 
otherwise  combustion  will  continue.  A  mathematical  model  was 
developed  to  provide  an  understanding  of  this  interaction 
and  its  effect  on  thermal  behavior. 

The  mathematical  model  is  formulated  as  follows.  Energy 
balances  on  the  carbon  and  convected  air  provide  heat 
transfer  equations  for  each.  The  heat  transfer  mechanisms 
incorporated  in  the  model  are:  (1)  conduction,  (2)  convection 
and  (3)  radiation  (where  applicable,  at  boundaries).  In 
addition,  nonvolatile  combustion  is  included  in  the  carbon 
heat  transfer  equation  as  a  fractional  order  heat  generation 
term  of  Arrhenius  type. 

The  conservation  of  species  law  is  applied  to  the  oxygen 
molecule  concentration  to  yield  a  third  equation.  The 
resulting  oxygen  mass  transfer  equation  includes  the 
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medium.  The  heat  transfer  mechanisms  included  were  conduc¬ 
tion,  convection  and  radiation.  One  problem  considered  was 
a  porous  mat  subject  to  Arrhenius  combustion  in  which  all 
properties  were  temperature  dependent.  The  capillary  serial 
(permeability)  model  introduced  by  Scheidegger  [Ref.  10]  was 
employed. 

The  current  investigation  considers  the  effects  of 
porosity  on  system  behavior  for  a  two-dimensional  (axi- 
symmetric)  model  of  combustion  in  a  porous  medium.  A  great 
part  of  the  following  analysis  may  be  found  in  Vatikiotis 
[Ref.  9].  The  pertinent  parts  are  repeated  here  for  conven¬ 
ience  to  the  reader.  Deviations  are  cited  and  are  for  the 
most  part  due  to  the  axisymmetric  two-dimensional  aspect  of 
the  present  model  vice  the  previous  one-dimensional  model. 

The  storage  scheme  employed  by  the  numerical  model  warrants 
the  reader's  attention.  The  savings  in  computer  storage 
realized  by  utilizing  the  method  of  Franke  and  Salinas 
[Ref.  29],  is  substantial  and  is  addressed  in  the  Numerical 
Section  of  this  work.  The  nonlinear  combustion  term  is 
treated  as  a  bilinear  spatial  operator  of  carbon  temperature 
and  oxygen  concentration,  in  contrast  to  the  Vatikiotis  treat 
ment  of  the  term  as  an  excitation  vector.  The  idea  behind 
the  present  treatment  was  to  capture  the  effects  of  the  non¬ 
linear  combustion  term  on  both  of  these  dependent  variables. 
It  was  felt  that  this  would  alleviate  some  numerical  diffi¬ 
culties.  A  numerical  model  is  formulated  and  results  are 
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The  boundary  configuration  for  this  work  was  a  two-dimensional 
region  with  a  line  of  vertical  symmetry,  enclosed  on  both 
sides  by  impermeable,  insulated  sidewalls,  and  heated  along 
half  of  the  base.  The  convective  flow  of  fluid  through 
porous  media  heated  from  below  has  applications  in  the  study 
of  behavior  of  geothermal  systems.  Hickox  studied  convection 
as  an  application  to  sub-seabed  disposal  of  nuclear  waste. 

The  problem  was  a  two-dimensional  transient  analysis  of  free 
convection  produced  by  a  concentrated  heat  source  (implanted 
container  of  waste  material)  in  the  subsurface  sedimentary 
layer  of  a  seabed.  The  porous  medium  was  assumed  to  be  rigid, 
homogeneous  and  isotropic.  Density  changes  of  the  fluid  were 
accounted  for  only  in  the  buoyancy  term  of  Darcy’s  law. 
Permeability,  viscosity  and  thermal  conductivity  were  assumed 
to  be  constant.  Chan  and  Banerjee  conducted  a  transient 
three-dimensional  analysis  of  natural  convection  in  porous 
media.  In  addition  to  convective  heat  transfer,  conduction 
was  also  considered  between  solid  spherical  particles  that 
comprise  the  porous  medium.  The  porous  medium  was  considered 
homogeneous  and  isotropic  in  its  physical  properties  including 
permeability  and  thermal  conductivity,  both  of  which  were 
assumed  to  be  constant  with  temperature.  Fluid  density  was 
considered  to  be  constant  except  in  the  buoyancy  term  of 
Darcy's  law. 

In  1982,  Vatikiotis  (Ref.  91  considered  a  transient  one¬ 
dimensional  heat  transfer  and  combustion  model  of  a  porous 
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Ignition  parameters  in  porous  solid  fuels  have  been 
analyzed  by  Kim  and  Chung  [Ref.  4].  They  investigated 
three  porous  solid  fuel  geometries  (a  semi-infinite  slab,  an 
infinitely  long  circular  cylinder  and  a  sphere)  with  con¬ 
stant  energy  and  gaseous  oxidant  fluxes  at  the  fuel  surface. 
Laplace  transformation  of  the  nondimensionalized  oxidant  mass 
equation  and  fuel  energy  equation  allowed  for  asymptotic 
solution  of  a  nondimensionalized  transient  temperature 
equation  which  is  valid  in  the  neighborhood  of  the  fuel 
surface.  Observations  included  shorter  ignition  times  for 
spheres  compared  to  slabs.  Times  for  ignition  increased  with 
fuel  size  and  approached  values  of  semi-infinite  slabs 
asymptotically.  Other  effects  of  size  and  geometry  of  porous 
solid  fuels  on  ignition  parameters  are  presented. 

Saatdjian  and  Caltagirone  [Ref.  5]  investigated  a  transient 
two-dimensional  combustion  model  with  natural  convection. 

A  porous  medium  undergoing  exothermic  combustion  was  saturated 
by  a  gas  and  bounded  by  two  impermeable  planes.  Permeability, 
fluid  viscosity,  thermal  conductivity  of  the  porous  medium, 
and  the  heat  of  reaction  were  assumed  to  be  constants. 

Horne  and  O'Sullivan  [Ref.  6],  Hickox  [Ref.  7]  and  Chan 
and  Banerjee  [Ref.  8],  investigated  the  natural  convection 
phenomenon  in  porous  media.  Exothermic  reactions  were  not 
addressed  in  these  investigations.  Horne  and  O'Sullivan 
considered  the  effects  of  variable  viscosity  on  the  stability 
of  a  porous  layer  in  a  transient  two-dimensional  problem. 


Investigations  by  Kordylewski  on  the  influence  of  aero¬ 
dynamics  on  the  critical  parameters  of  thermal  ignition 
[Ref.  1],  examined  homogeneous  combustion  in  a  two-dimensional 
cylindrical  reactor.  The  transient  analysis  involved 
Arrhenius  combustion  of  a  solid  fuel.  Heat  transfer  mechanisms 
considered  were  convection  and  conduction.  Because  the  inves¬ 
tigation  dealt  with  thermal  ignition  theory,  reactant  con¬ 
sumption  was  omitted.  The  flow  field  was  assumed  to  be 
steady  prior  to  ignition  and  constant  fluid  properties  were 
assumed. 

Safety  dictates  the  assessment  of  the  structural  integrity 
of  a  building  after  a  severe  fire.  In  order  to  predict 
stresses  due  to  fire  in  buildings,  the  thermal  response  must 
be  known.  Sahota  and  Pagni  [Ref.  2] ,  formulated  a  transient 
solution  for  two-phase,  two-component  flow  in  one-dimensional 
porous  concrete  structures.  The  mechanisms  considered  in¬ 
cluded:  heat  conduction,  molecular  diffusion  of  gaseous 

components,  and  pressure-driven  convective  flow  subject  to 
Darcy's  law. 

A  sudden  reduction  in  feed  temperature  in  a  packed-bed 
reactor  leads  to  the  transient  temperature  rise  known  as 
"wrong-way  behavior."  This  behavior  was  investigated  by 
Mehta,  Sams,  and  Luss  [Ref.  3].  The  work  identifies  the 
important  rate  processes  and-  parameters  which  cause  the 
behavior,  and  generates  an  expression  for  predicting  the 
magnitude  of  the  maximum  transient  temperature. 
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I. 


INTRODUCTION 


A.  PRIOR  INVESTIGATIONS 

The  problem  of  characterizing  physical  systems  involving 
combustion  and  quantifying  the  accompanying  thermal  response 
has  received  attention  in  recent  years.  The  level  of  diffi¬ 
culty  encountered  in  a  problem  of  this  type  precluded 
analytical  as  well  as  numerical  solutions.  The  difficulty 
of  the  problem  lies  in  the  number  of  disciplines  that  encom¬ 
pass  it.  The  problem  involves  the  kinetics  of  combustion, 
heat  and  mass  transfer  mechanisms  and  fluid  flow.  Numerical 
solutions  with  increased  efficiency  in  computation  are  now 
possible  with  high  speed  computers. 

The  combustion  problem  has  applications  in  the  areas  of 
forest  fire  control,  energy  conservation,  underground 
(nuclear)  fuel  storage  and  waste  disposal,  and  natural  gas 
fire  control.  Major  contributions  may  be  made  in  the  area 
of  energy  production  and  conservation  by  the  analysis  of 
heat  generation  and  ignition  characteristics  of  combustible 
materials . 

The  characterization  of  combustion  and  heat  transfer  in 
porous  media  has  been  of  particular  interest.  A  brief  survey 
of  some  of  the  investigations  in  the  area  of  combustion  and 
heat  transfer  in  porous  media  follows.  The  intent  is  to 
acquaint  the  reader  with  work  in  the  field  that  offered 
insight  into  the  formulation  of  this  investigation. 
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K 

5  --  Kronecker  delta  function 

£  --  Thermal  emissivity 

6  --  Solution  coefficient 

A  —  Field  operator 

U  --  Dynamic  viscosity 

S  --  Local  element  coordinate 

j  --  Mass  density 

j  --  Stef an-Boltzman  constant 

t  --  Tortuosity,  stress 

p  --  Oxygen  concentration 

V  —  Particle  shape  factor,  approximate  solution 


CO2  —  Carbon  dioxide 
e  --  Effective 


fm  --  Film 


g  --  Heat  generation 
i  --  At  the  current  time  or  step 
ig  —  Ignition 

o  —  Cylinder  dimension  (i.e.,  rQ  is  cylinder  radius, 
zq  is  cylinder  length) 

0^  --  Oxygen 

p  --  At  constant  pressure 
r  --  Radiation 
s  —  Solid 
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constitutive  property  (i.e.,  a  property  specified  by  a 
constitutive  relation,  Darcy's  law).  The  specific  permea¬ 
bility  is  proportional  to  the  filter  velocity  induced  by  a 
unit  pressure  gradient.  Darcy's  law  is  not  derived  from 
first  principles  but  through  exhaustive  experimental 
analyses . 

The  Dupuit-Forcheimer  assumption,  addressed  in  Carman 
[Ref.  11],  relates  the  loca±  pore  velocity  to  the  filter 
velocity  by, 

Q  =  p  V  (3.7) 

The  hypothesis  is  that  the  local  pore  velocity  is  greater 
than  the  filter  velocity.  The  actual  velocity  in  a  single 
pore  is  a  function  of  the  fluid  element's  location  within 
the  pore.  The  Dupuit-Forcheimer  assumption  defines  an 
"average"  velocity  within  the  pore. 

The  continuity  equation  for  two-dimensional  flow  in 
porous  media  with  nonconstant  porosity  distributions  is, 


°(PP  J 

— — +  p  p  Div  V  =  0  (3.8) 

Dt  a 


Invoking  the  Dupuit-Forcheimer  assumption  and  Darcy's  law, 
the  continuity  equation  becomes, 


D  (pp 


n .  , -m  ,  3P 

p  p  Div  — (-r— 
^  a  pp  3r 


,3P 

(3l  “  p  pagU) 


=•  0 


(3.9) 
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Equation  3.9  'neglecting  body  forces)  is  one  of  four  field  equations 
cast  into  a  finite  element  formulation  later  in  this  work.  From  the 
pressure  distribution,  the  pore  velocity  distribution  is 
obtained  by  invoking  Darcy's  law  and  the  Dupuit-Forcheimer 
assumption.  The  derivation  of  Equation  3.9  is  presented  in 
Appendix  A. 

C.  SEMENOV  MODEL  OF  COMBUSTION 

The  model  of  Semenov  described  in  Frank-Kamenetskii 
[Ref.  13]  and  Vulis  [Ref.  14],  is  the  combustion  model 
employed  herein.  The  basis  of  the  model  is  the  relation  of 
reaction  rate  to  temperature  and  the  interaction  of  heat 
generation  and  heat  transfer.  The  reaction  rate  expression 
Rc ,  is  the  Arrhenius  expression  for  a  simple  n-th  order 
reaction , 


R  =  A  ^nexp( — p— ) 

R  T 
u  c 


(3.10) 


where  A  is  the  time  constant  of  the  chemical  reaction,  E  is 
the  activation  energy,  R^  is  the  Universal  G^s  Constant, 

Tc  is  the  absolute  carbon  temperature,  and  <p  is  the  oxygen 
concentration.  In  a  simple  reaction,  the  reaction  rate 
depends  on  the  concentrations  of  reactants  and  not  on  the 
products.  The  heat  generated  by  the  reaction  is  obtained 
by  multiplying  Expression  3.10  by  the  heat  (enthalpy)  of 
combustion.  As  explained  in  [Ref.  13],  a  theory  of  combustion 
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can  be  constructed  only  if  certain  reasonable  assumptions 
are  made.  For  example,  in  order  to  regard  the  initial 
states  as  stationary,  one  must  neglect  the  reaction  rates 
at  these  states.  The  empirical  law  of  chemical  kinetics 
embodied  in  the  Arrhenius  expression  tells  us  that  the  reac¬ 
tion  rate  never  goes  to  zero  but  falls  off  exponentially 
with  a  decrease  in  temperature.  Neglect  of  the  reaction 
rate  at  the  initial  states  is  necessary  to  achieve  initial 
equilibrium.  In  a  finite  range  of  temperatures  neve  me 
initial  states,  neglect  of  reaction  rates  is  also  justifiable 

[Ref.  13].  At  room  temperature  the  reaction  rates  are  on 

2 

the  order  of  l.E-15  lbm-carbon/f t  -hr  or  smaller,  vice 

2 

l.E+3  lbm-carbon/f t  -hr  at  1500  °F. 

D.  ARRHENIUS  LAW  OF  REACTION  RATES 

In  1934,  Parker  and  Hottel  [Ref.  15]  proposed  an 
Arrhenius  expression  for  the  reaction  rate  of  carbon  reacting 
in  air.  The  expression  assumed  a  simple  first  order  reaction 
for  the  combustion  of  carbon  in  air.  Frank-Kamenetski i 
[Ref.  13]  has  shown  that  the  Parker  and  Hottel  data  is  better 
correlated  by  a  fractional  order  reaction,  n  between  1/3 
and  2/3.  A  reaction  order  of  1/2  yields, 

r  1/2 

R  =  2.065  x  10  (Rn  t>)  /  exp  (-57, 2  40/R  -T  )  (3.11) 

c  u  c 

2 

In  Expression  3.11,  Rc  is  in  units  of  lbm-carbon/f t  -hr, 

Rq  is  the  gas  constant  for  oxygen  (48.29  ft-lbf/lbm-R) , 
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S  is  in  lbm/ft^,  is  1.986  Btu/1 bmole-R,  and  Tc  is  in 
Rankine.  In  order  to  determine  the  rate  of  heat  generation 
and  the  rate  of  oxygen  consumption,  the  chemical  reaction 
for  the  combustion  process  must  be  considered.  Although  many 
complex  chains  in  chemical  kinetics  generally  describe 
combustion,  a  simple  two-product  analysis  is  employed  in  this 
formulation.  For  nonvolatile  combustion  of  carbon  and 
oxygen,  two  reactions  that  describe  the  process  are, 


C 


CO 


(3.12) 


C  +  O2  CO2 


(3.13) 


where  0  denotes  oxygen;  and  C,  CO  and  C02  denote  carbon, 
carbon  monoxide  and  carbon  dioxide,  respectively.  The  ratio 
of  the  mass  rates  of  carbon  monoxide  to  carbon  dioxide 
produced  increases  with  increasing  temperature.  Arthur 
[Ref.  16]  presents  an  expression  for  the  rate  ratio  as  a 
function  of  temperature  (in  Kelvin). 


~~  =  2500  exp(-6240/Tj 


(3.14) 


The  expression  is  valid  for  temperatures  between  790  and 
1170  K  (1310-2110  degrees  Fahrenheit).  As  a  result  of  this 
temperature  dependency,  the  s tochiometr ic  ratio  and  the  heat 
of  reaction  are  functions  of  temperature.  Defining  the 
fraction  of  carbon  monoxide  being  produced  by » 
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(3.15) 


CO 


,  cu 
lCO„ 


)/ (1  + 


(— )  ) 
vCO  '  ' 


and  the  fraction  of  carbon  dioxide  as, 


CO. 


=  1/(1  +  ( 


CO 

C0- 


) ) 


the  heat  of  combustion  is  then  expressed  as, 


(3.16) 


OH 


R 


Fco  _Hco  +  fco2  _hco2 


(3.17) 


Values  for  the  heats  of  combustion,  uHCQ  and  ,  as 

functions  of  temperature  may  be  obtained  from  JANAF  (Joint 
Army,  Navy,  and  Air  Force)  tables  [Ref.  17].  For  the  range 
of  temperatures  being  investigated  (80-2000  degrees  Fahren¬ 
heit)  ,  the  heats  of  combustion  are  3966.3  Btu/lbm  for  FCQ 

and  14,121.  Btu/lbm  for  F  .  Frank-Kamenetskii  [Ref.  13] 

^  ^  2 

points  out  that  over  narrow  ranges  of  temperature  it  is 
permissible  to  use  an  approximate  expression  which  correctly 
describes  the  reaction  rate.  The  stochiometr ic  ratio  (fuel 
to  oxygen)  of  the  overall  reaction  is. 


■R 


fCO  fCO. 


/(f 


C0o  CO 


f  F 
CO  CO. 


(3.18) 


where  f is  the  stochiometric  ratio  for  the  reaction  3.12 
and  f  is  the  stochiometric  ratio  for  the  reaction  3.13. 
The  rate  of  heat  generation,  R  , 


and  the  rate  of  oxygen 


consumption  may  now  be  expressed  by, 

(3. 19] 

(3.20) 

Parker's  and  Hottel's  work  [Ref.  15]  and  Arthur's  work 
[Ref.  16]  were  conducted  with  specific  types  of  carbon. 

Tables  and  references  exist  in  Smoot  and  Pratt  [Ref.  18]  and 
Frank-Kamenetskii  [Ref.  13]  for  rate  expressions  using  other 
types  of  carbon.  The  present  model  allows  for  any  simple 
fractional  order  rate  expression  of  Arrhenius  type  to 
account  for  carbon  consumption  with  CO  and  CO^  byproducts. 

For  this  work  the  Parker  and  Hottel  rate  expression  as  modi¬ 
fied  by  Frank-Kamenetskii  (n  =  1/2)  is  used. 

Particle  diameter  decreases  with  progressive  combustion. 
The  rate  of  decrease  depends  on  the  amount  of  carbon  con¬ 
sumed  at  a  point  over  time.  Observations  have  shown  that 
the  effect  of  decreasing  particle  diameter  is  significant 
when  the  reaction  is  concentrated  in  a  small  region  of  the 
porous  medium.  To  account  for  this,  an  expression  for  the 
time  rate  of  diameter  change  as  a  function  of  reaction  rate 
is  derived  (Appendix  B) .  For  spherical  particles  the  equation 
is, 


O 

d 


- 2 R  Z  D3/(tto  d2; 
c  c 


(3.21) 


31 


where  .;  is  the  bulk  mass  density  of  carbon.  The  diameter 
and  the  reaction  rate  are  functions  of  time  and  space. 

E.  CARBON  HEAT  TRANSFER  EQUATION  FOR  POROUS  MEDIA 

The  heat  transfer  equations  of  the  present  investigation  incor¬ 
porate:  (1)  radiation  (at  boundaries  where  applicable),  (2) internal 

convection,  (3)  conduction,  (4)  internal  combustion,  (5)  temperature 
dependency  of  properties,  and  (6)  compressibility  effects 
of  air  into  a  two-dimensional  (cylindrical  coordinate) 
formulation.  The  carbon  energy  conservation  equation  is, 

IT 

(l-p)ccCc  =  7-  (1-p)  (ke)  (7Tc)  -hZ(Tc-Ta)  +  RgZ  (3.22) 


The  derivation  of  Equation  3.22  is  presented  in  Appendix  A. 
The  effective  conductivity,  k  ,  of  the  porous  solid  was 
proposed  by  Russel  [Ref.  24], 


k 

e 


,2/3 


+  -  P 


.2/3, 


,2/3 


-  p'  +  ir(i  -  ? 

c 


,2/3 


+  P') 


(3.23) 


where  k  and  k  are  bulk  thermal  conductivities  of  carbon 
c  a 

and  air  respectively,  and  p'  =  1-p.  Russel's  expression, 
which  is  based  on  an  electrical  analogy,  is  valid  for  the 
full  range  of  porosities,  0.0  to  1.0. 

Because  of  the  difficulties  encountered  in  a  radiative 
analogy  to  Fourier's  law  of  conduction,  the  particle  to 
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particle  radiative  exchange  in  the  porous  medium  is  omitted. 

The  difficulties  are:  (1)  the  geometry  does  not  easily  allow 
one  to  derive  an  expression  for  a  "sink"  temperature  to  use  in 
a  linearized  approximation  of  the  Stefan-Boltzman  equation,  and 
(2)  the  multi-wavelength  characteristics  of  the  radiation 
phenomenon  are  not  easily  incorporated. 

F.  HEAT  TRANSFER  EQUATION  FOR  AIR  IN  POROUS  MEDIA 

The  internal  convective  heat  transfer  coefficient  of 
Yoshida,  Ramaswami,  and  Hougen  [Ref.  33],  is  given  by, 

h  =  0.91  Re,_0-51[^C  C(C  y/k  )~2/3]f  (3.24) 

a  a  a  m 

where  :i>  is  equal  to  1  for  spherical  particles,  G  is  a  pseudo 

mass  velocity  given  by  p  o  s  and  C  is  the  specific  heat 

a  a 

of  air  at  constant  pressure.  The  f  subscript  indicates 
properties  are  to  be  evaluated  at  film  temperature.  Re'  is 
a  pseudo  Reynolds  number  defined  by, 

Re'  =  G/z  y  \\>  (3.25) 

The  air  properties,  as  well  as  the  internal  convection  coeffi¬ 
cient,  h,  are  temperature  dependent.  The  reaction  rate  in 
Equation  3.22  is  given  by  Expression  3.19. 

An  energy  balance  on  the  air  within  the  porous  medium 
obtains  the  second  heat  transfer  equation. 
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(3.26) 


po  C 


a  a  jt 


:•  (p  k  7T  )  +  hZ  (T  -T  ) 
a  a  c  a 


-  P  Paca(V'7)Ta  ~  (V‘7)P  p 


The  density  of  air,  c  ,  is  approximated  by  the  ideal  gas  law. 
Pressure  terms  are  due  to  the  compressibility  of  air.  The 
derivation  of  Equation  3.26  is  presented  in  Appendix  A.  All 
properties  in  Equation  3.26  are  temperature  dependent.  The 
properties  of  standard  air  were  used  in  the  model.  Vatikiotis 
[Ref.  9]  points  out  that  tolerable  differences  (average  of 
7%  difference)  between  standard  air  properties  and  properties 
accounting  for  the  presence  of  byproducts  CO  and  CO^ /  is 
acceptable.  Increased  accuracy  would  introduce  an  additional 
mass  balance  equation  for  either  CO  or  C02-  Polynomial  ex¬ 
pressions  used  to  calculate  the  thermophysical  properties  of 
air  are  those  presented  by  Vatikiotis  [Ref.  9].  The  expressions 
are  presented  in  Appendix  B. 

G.  OXYGEN  DIFFUSION  EQUATION  FOR  POROUS  MEDIA 

The  fourth  field  equation  necessary  to  complete  the 
system  of  equations  is  provided  by  an  oxygen  species  conser¬ 
vation  requirement.  Transport  mechanisms  included  in  the 
model  are  molecular  diffusion  (Fick's  law),  convective  mass 
flow  and  oxygen  consumption  due  to  combustion.  Pressure  and 
temperature  induced  concentration  gradients  are  considered 
negligible.  Vatikiotis  [Ref.  9],  provides  examples  for  which 
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diffusion  arising  from  pressure  and  temperature  gradients 
is  important.  The  oxygen  mass  balance  equation  is, 


P 


at 

j  t 


h  (p  ;))  -  v  •  (p  ?  V) 


(3.27) 


The  derivation  of  Equation  3.27  is  presented  in  Appendix  A. 
The  effective  diffusivity,  V  ,  proposed  by  Denbigh  and 
Turner  [Ref.  31]  for  a  porous  medium  is 


De  =  V/x 


(3.28) 


Expression  3.28  accounts  for  the  tortuosity  encountered  by 
the  oxygen  molecules  as  they  flow  through  the  porous  medium. 
The  semi-empirical  expression  proopsed  by  Gilliland  [Ref.  34] 
is  used  to  obtain  the  diffusivity  of  oxygen  into  air.  The 
expression  is, 


V  =  435.7  T 


3/2  ( m"1  +M"1)  1/2/[P(Vy3  +Vq/3)2J 

a  cl  ^2 


(3.29) 


2 

where  V  is  in  units  of  cm  /s,  P  is  the  total  pressure  in  Pa, 

V  and  V  are  the  molecular  volumes  of  air  and  oxygen, 

^  u  2 

respectively.  M  and  M  are  the  molecular  weights  of  air 

a  o2 

and  oxygen.  The  values  of  V  and  V  may  be  obtained  in 

<31  U  ^ 

3Z 

Holman  [Ref.  35]  as  29.9  and  7.4  cm  ,  respectively.  The 
oxygen  consumption  term  in  Equation  3.27  is  given  by  Expression 
3.20. 
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H.  BOUNDARY  CONDITIONS 

Boundary  conditions  employed  were  as  follows  for  the 
carbon , 
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5r  -L«o 

1  r 
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=  0 


(3.30) 
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=  -3, 
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(3.31) 


where  qg  is  the  starting  heat  flux. 
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(1-D  problems) 
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(3.35a) 
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For  pressure, 
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Equations  3.35,  3.36,  3.43,  and  3.44  correspond  to  convective 

flux  conditions  on  air  temperature  and  concentration. 

At  the  air  inlet  (—  =  0.0),  Dinchlet  conditions,  Equations 

L  o 

3.35a  and  3.43a  may  be  imposed  on  the  air  temperature  and 
oxygen  concentration.  The  idea  behind  this  treatment  is  that 
in  the  presence  of  a  semi-infinite  medium  (ambient  air) ,  the 
air  and  0^  concentration  may  be  considered,  to  a  first 
approximation,  very  near  ambient  conditions.  Equations  3.'0, 
3.34,  3.28  and  3.42  correspond  to  symmetry  conditions  at 
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=  0.0.  Equations  3.33,  3.37,  3.41  and  3.45  correspond 
“o 

to  impermeable  boundaries  and  insulated  boundaries  at  —  =  1.0 

ro 

and  are  used  in  conjunction  with  one-dimensional  problems. 

The  second  part  of  Equations  3.33  and  3.37  are  the  boundary 
conditions  that  represent  a  relaxation  of  the  radial  insulation 
of  carbon  and  air  temperatures  which  allow  the  system  to 
become  a  two-dimensional  heat  transfer  problem.  It  must  be 
pointed  cut  however  that  the  heat  transfer  coefficients  of 
Equations  3.23  and  3.37  are  not  easily  obtained.  This  work 
could  net  proceed  if  the  following  assumptions  were  not  made. 

If  was  assumed  that  the  air  and  carbon  had  near  equal  profiles 

at  ~  =  1.J  so  that  the  same  heat  transfer  coefficient  would  apply 

o 

to  both  variables.  If  Equation  3.35  is  used  as  a  boundary 
condition,  the  air  temperature  follows  closely  the  carbon  tempera¬ 
ture  at  — -  =  1.0.  Furthermore  it  was  assumed  for  simplicity  cf 

“  Q 

calculation  that  the  heat  transfer  coefficient  was  constant. 

The  difficulty  is  as  follows.  In  order  to  determine  the 

neat  transfer  coefficient  for  the  impressed  temperature 

r  l-.-s  at  ~  =  1.0,  the  temperatures  themselves  must  be 
"  o 

•ir.  wr. .  !  correct  procedure  would  involve  an  initial  esti¬ 

mate  ::  the  profile,  calculation  of  the  heat  transfer  coeffi- 
i  :-*nt,  solution  of  the  problem  for  one  integration  and 
veri float icn  of  the  assumed  and  calculated  temperature 
r  ;f:ies.  Upon  suitable  verification  of  the  profiles,  the 
same  procedure  is  performed  for  the  next  integration.  Other 
problems  arise  for  example  if  the  radial  heat  transfer  at 
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c . 

1 


(4.14) 


4  sin  ■ 

2  3.  +  s in ( 2  2 • ) 

l  l 

The  constants  2 ^  are  the  roots  of  the  transcendental  alge¬ 
braic  equation, 

2.  tan  3  .  =  B.  =  h  L/k  (4.15) 

11  1  o 

Thus  the  solutions  have  the  Biot  number  of  the  slab  as  a 
parameter.  The  roots  of  Equation  4.15  are  tabulated  in 
Appendix  A  to  [Ref.  25] . 

3.  The  Cvlinder 

■*  -  - 

Consider  the  cylinder  of  Figure  4.1(b).  The  suddenly 
immersed  cylinder  satisfies  the  equation, 


3T 
j  t 


1 

r 


(4.16) 


subject  to  an  initial  condition, 
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5.E-1  hr.,  and  iterated  about  the  given  initial  fields. 

This  sufficiently  demonstrated  the  first  of  several  tests 
requisite  to  assure  numerical  code  compatibility  with  the 
integration  scheme. 

2 .  The  Finite  Slab 

Consider  the  slab  of  thickness  2L  in  Figure  4.1(a). 
One  seeks  the  solution  of  the  one-dimensional  conduction 
equation , 


•? 

•  T  -■>  “T 


subject  to  an  initial  condition, 


(4.10) 


T  ( x ,  0  )  -  T  i 


(4.11) 


and  uniform  convection  conditions  at  both  surfaces:  At 
x  =  *.  L : 

-k  4—  =  ±  h (T  -  T  )  (4.12) 

jX  o 

The  exact  solution  is  given  as  a  Fourier  series  in  [Ref.  23], 


T  -  T 
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T  ■  -  T 
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C  .  e 
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1  L 


(4.13) 


where , 
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2  2  -22  nodal  points.  The  solution  changed  5%  from  a  12  -12 
to  a  22  ■  22  nodal  point  model.  Next,  a  17  <  17  point  non- 
uniform  grid  model  was  investigated.  It  was  found  that 
the  solution  changed  approximately  2%  from  the  22  <  22  nodal 
point  model.  At  this  point,  it  was  determined  that  the 
12  ■  12  non-uniform  nodal  point  afforded  the  desired  balance 
between  cost,  computational  effort  and  computer-run  time. 


grid  size 

kmax 

clock 

runtime  (min . ) 

(n  *  n) 

12 

6  ,  398 

70 

17 

13,223 

100 

22 

22,498 

110 

MODEL  VALIDATION  TESTS 

In  order  to  validate  the  capabilities 

and  accuracy  of 

the  numerical  code,  several  tests  were  conducted. 

1 .  The  Steady  State  Problem 

First,  it  was  necessary  to  ensure  the  model's  ability 
to  recognize  a  steady  state  condition.  This  was  the  simplest 
of  all  tests  (and  first  to  be)  performed.  Initial  (ambient) 
conditions  were  imposed  uniformly  throughout  the  system  on 
the  four  fields:  carbon  temperature  and  air  temperature 
(80  degrees  F.),  pressure  (2,116.8  psf)  and  oxygen  concen¬ 
tration  (0.0172  lbm/ft^).  An  initial  integration  time  step 
equal  to  the  minimum  time  step  allowable  (user  input)  l.E-6  hr. 
was  selected.  In  less  than  five  integrations  the  algorithm 
switched  to  the  maximum  allowable  time  step  (user  input) 
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superior  to  uniform  grids.  In  fact,  for  many  problems, 
computer  storage  limitations  require  that  a  non-uniform  grid 
be  used  to  obtain  reasonably  accurate  solutions.  N’on- 
uniform  grids  allow  the  investigator  to  take  advantage  of 
the  knowledge  of  areas  of  severe  combustion  induced  gradients. 
"Stacking"  elements  in  these  areas  assists  the  algorithm  in 
producing  more  accurate  and  numerically  stable  results.  For 
similar  degrees  of  accuracy  with  uniform  grids,  it  is  esti¬ 
mated  that  grids  on  the  order  of  1000  *  1000  nodal  points 
(and  larger)  would  have  been  required.  Nondimensionalized 
length  discretizations  on  the  order  of  .001-. 005  and  smaller 
near  the  air  inlet  boundary  permit  excitations  on  the  system 
boundaries  to  be  propagated  accurately  through  the  medium. 
Until  discretizations  on  the  order  of  non-dimensionalized 
lengths  of  0.001  were  employed,  numerical  instability  was 
encountered  in  the  carbon  temperature  field.  This  instability 
was  manifested  by  severe  overshoot,  leading  to  negative 
temperatures  in  close  proximity  to  severe  thermal  gradients. 

For  non-uniform  grids  it  was  found  that  grids  below  10 
nodal  points  in  each  direction  were  not  adequate.  The 
temperatures  obtained  from  these  models  were  lower  due  to  the 
relatively  low  degree  of  grid  refinement  or  discretization 
possible  with  such  few  points.  Grid  refinement  in  this  type 
of  problem  is  essential  to  capture  the  high  activity  (large 
gradients)  that  is  inherent  in  the  combustion  phenomenon. 

A  convergence  study  was  done  on  several  nodal  points  at 
various  times  for  non-uniform  grids  at  12  <12,  17  <17,  and 
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combustion  occurring  within  the  element.  An  element  loop  is 
performed  to  account  for  combustion  terms  in  an  elementally 
averaged  sense  over  the  entire  system  spatial  domain. 

B.  OPTIMAL  COMPACT  STORAGE  (OCS) 

The  optimal  compact  storage  scheme,  presented  by  Franke 

and  Salinas  [Ref.  29],  employed  in  the  solution  of  the 
combustion  problem,  allowed  the  computational  effort  to 
proceed  with  substantial  savings  in  computer  storage, 
computer-run  time,  and  ultimately  dollar  cost  per  run.  In 
the  optimal  compact  storage  scheme,  only  the  non-zero  coeffi¬ 
cients  are  stored  in  a  vector  array.  This  results  in  a 
very  significant  reduction  of  the  storage  area  compared  to 
banded  storage.  One  might  encounter  problems  on  the  order 
of  1000  <1000  DOF.  (In  this  problem  there  are  four  degrees 
of  freedom  at  each  nodal  point.  A  17  x 17  grid  requires  a 
1156  < 1156  matrix  if  full  storage  is  employed.)  For  banded 
storage  the  bandwidth  might  be  200.  The  storage  ratio  for 
bandwidth  to  full  storage  would  be  0.2.  Using  OCS,  a  con¬ 
servative  estimate  of  the  average  number  of  equation  entries 
in  this  model  is  20.  The  ratio  of  OCS  to  banded  storage 
is  0.1.  Thus  in  this  example,  the  savings  realized  by  OCS 
vice  full  storage  is  98%  ! 

C.  GRID  CONVERGENCE 

A  convergence  study  was  done  on  the  response  field  for 
several  grids.  It  was  found  that  non-uniform  grids  were  far 
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3.  Bilinear  ODerator  Treatment  of  0_  and  T 

- * - 2 - c 

The  bilinear  operator  treatment  of  the  reaction  rate 


is  the  present  method  of  treatment.  The  reaction  rate  term 


is  rearranged  as  follows, 


A  :n  exp (-E/RT  )  *t  . 

R  = - 1 - S-:  "-1  •  [T  ,1  (4.6) 

C  T  c 

c 


Letting 


T 


c 


irli 


i  =  1,3 


(4.7) 


and 


_T\ 

’ir4i  ' 


i  =  1,3 


(4.8) 


and  invoking  natural  coordinates  [Ref.  28],  an  elementally 
averaged  contribution  results  in  the  following  area  integral, 

C  =  C  /  /  C  tT9.  dA  (4.9) 

Ae . 

The  double  subscript  permutation  of  carbon  temperature  and 
C>2  concentration  results  in  a  3x9  elemental  matrix  that  is 
distributed  into  the  system  equations.  In  this  scheme  at 
each  nodal  point  within  an  element,  carbon  temperature  and 
concentration  equations  receive  nine  terms  arising  from 
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In  the  carbon  equation,  Expression  4.1  appears  as, 


R  =  !H_  •  R_ 

g  R  c 


In  the  oxygen  concentration  equation,  Expression  4.1  appears 


R02  f  R  '  Rc 


The  term  that  is  to  be  evaluated  at  the  last  time  step  and 
to  be  held  constant  over  the  next  time  step  is, 


Rc  =  {A  4>  nexp(-E/RTc)  } 


where  the  superscript  *tn_^  indicates  evaluation  occurring 
at  the  previous  time  step. 

2.  Linear  Operator  Treatment  of  0-.  Concentration 


In  order  to  realize  an  improvement  over  the  first 
treatment  one  may  retain  a  portion  of  the  reaction  expression 
as  an  operator  by  making  the  following  rearrangement: 


Rc  =  {A  bn_1 exp(-E/RTc) }  n_1 


(4.5) 


where  [5]  denotes  a  spatial  operator  treatment  of  the 
response  variable. 
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IV.  NUMERICAL  CONSIDERATIONS 


A.  TREATMENT  OF  NONLINEAR  COMBUSTION  TERM 

The  treatment  of  the  nonlinear  combustion  term  is  now 
discussed  in  brief  detail  highlighting  the  advantages  of 
each  treatment.  A  more  comprehensive  discussion  is  presented 
in  Appendix  C.  The  interested  reader  is  encouraged  to 
review  the  detailed  analysis.  There  are  several  ways  to 
treat  the  Arrhenius  reaction  rate  expression. 

1 .  Excitation  (Force)  Term  Treatment 

Perhaps  the  simplest  method  for  incorporating  the 
combustion  term  is  to  evaluate  it  at  each  time  step  and  use 
the  evaluated  value  (held  constant)  as  an  estimate  of  the 
mean  value  during  the  next  integration.  It  is  realized  that 
in  fact,  the  value  is  not  constant  in  real  time.  This  treat¬ 
ment,  however,  provides  a  means  of  incorporating  the  term 
into  the  system  of  equations  as  a  first  approximation  with 
relatively  little  computational  effort.  Averaging  techniques 
exist  for  improving  upon  this  method.  In  this  scheme  the 
excitation  vector  is  modified  at  each  nodal  point  in  the 
carbon  temperature  and  oxygen  concentration  equations  to 
include  the  combustion  terms.  The  exponential  reaction  rate 
term  that  appears  in  both  the  porous  solid  and  oxygen  diffu¬ 
sion  equation  is, 


Rc  =  A  fn exp(-E/RTc) 


(4.1) 
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lectrical  Analogy  of  Thermal  Circuit 
t  r/rQ  *1.0 


5 


of  Porous  Medium 


Expression  3.46  is  incorporated  into  the  numerical  formula¬ 
tion  as  Neumann  boundary  condition  for  carbon  temperature. 
The  initial  heat  flux  may  be  turned  off  at  any  specified 
time.  The  boundary  condition  of  Expression  3.46  then 
switches  to  a  Cauchy  (mixed)  boundary  condition  to  account 
for  radiative  heat  transfer  from  the  carbon  particles  at  the 
boundary  to  the  ambient  air.  Convective  heat  transfer  from 
the  carbon  particles  to  the  air  is  treated  through  the 
internal  heat  transfer  coefficient.  The  above  procedures 
for  treating  problem  initiation  obviate  the  need  of  trying 
to  specify  initial  conditions  which  may  in  general  be 
arbitrary  for  each  problem. 
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analog  sense  (Figure  3.3)  and  thus  affects  the  temperature 
profile  and  ultimately  the  heat  transfer  coefficients.  Of 
necessity  then,  the  material's  inherent  ability  to  absorb 
energy  is  denied.  No  attempt  is  made  here  to  conjecture  as 
to  how  physically  one  may  impose  the  radial  boundary  conditions 
implemented  below.  If  it  were  possible  to  rigorously  apply 
the  heat  transfer  boundary  condition  described  above,  the 
numerical  model  would  be  capable  of  obtaining  a  solution. 

The  model  was  tested  at  various  values  of  constant  heat  trans¬ 
fer  coefficient:  0,  1,  5,  50,  100,  500,  to  determine  the 
effect  of  the  heat  transfer  coefficient  on  the  solution.  It 
was  found  that  until  the  heat  transfer  coefficient  gets 
significantly  large  compared  to  the  start  flux  applied  to 

the  carbon  at  the  air  inlet  boundary  (in  this  case  above 
2 

100  Btu/ft  -hr  compared  to  1500  flux  applied  to  carbon) ,  the 
two-dimensional  effects  on  the  temperature  profile  were 

r* 

localized  in  the  —  =  0.85  to  1.0  region.  The  temperature 

o 

IT  ... 

profiles  for  —  less  than  0.85  were  essentially  one-dimensional 
o 

(constant  value  independent  of  r) . 

I.  INITIAL  CONDITIONS 

The  model  is  developed  so  each  problem  begins  at  a 
uniform  initial  temperature.  A  heat  flux  applied  to  the 
carbon  along  a  boundary  provides  a  means  of  bringing  the 
system  to  a  temperature  level  where  reaction  terms  are 
sufficiently  high  to  generate  combustion.  The  heat  flux  is 
treated  as  follows, 
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—  =  1.0  is  of  natural  convective  type.  At  the  exit  —  =  1.0, 
r  z 

o  o 

functional  and  derivative  continuity  must  be  demonstrated 
for  pressure,  oxygen  concentration  and  air  temperature  at 
the  vertical  boundary  separating  the  stream  of  exiting  air 
and  the  convective  boundary  layer.  This  is  known  as  a  conju¬ 
gate  problem.  The  above  method  is  a  realistic  method  for 
-_he  treatment  of  this  problem  but  admittedly  it  is  beyond  the 
scope  of  this  work.  The  numerical  code  has  provisions  for 
incorporating  an  average  of  isoflux  and  isothermal  natural 
convective  Nusselt  number  formulations  (Churchill  correlations) 
for  vertical  cylinders  obtained  from  [Ref.  36].  Although  the 
subroutine  has  been  written  there  has  been  no  attempt  to 
date  to  implement  it.  It  is  known  that  physically  heat 
transfer  coefficients  lie  between  those  generated  from  iso¬ 
thermal  and  isoflux  considerations.  However,  simply  having 
empirical  correlations  for  evaluating  the  heat  transfer 
coefficients  does  not  eliminate  the  need  of  the  iterative 
scheme  described  above.  Additional  considerations  must  be 
addressed  in  this  context.  In  order  to  effect  an  impermeable 

vertical  boundary  at  —  =  1.0,  the  system  must  be  enveloped 

o 

by  a  cylindrical  container  fabricated  of  some  type  of  material. 
This  material  has  an  inherent  potential  to  absorb  internal 
energy  that  transits  from  the  initial  system  to  the  boundary 
where  some  flowing  medium  is  able  to  convect  energy  away 
(Figure  3.2).  The  material's  ability  to  absorb  energy  adds 
a  thermal  resistance  to  the  heat  transmission  in  an  electical 
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(4.19) 


2iJl^i)/Jo(3i)  =  Bi  =  hoVk  (4'21 

The  functions  J  and  J.  are  the  Bessel  functions  of  the  first 
o  1 

kind  whose  numerical  values  are  tabulated  in  most  advanced 

engineering  mathematics  texts.  For  a  limited  range  of 

interest,  numerical  values  of  J  and  J,  are  tabulated  in 

o  1 

Appendix  A  of  [Ref.  25].  The  roots  of  Equation  4.21  are 
tabulated  in  [Ref.  23]. 

4 .  Multidimensional  Solutions  by  the  Product  Method 

The  classic  solutions  for  semi-infinite-  and  finite¬ 
thickness  bodies  (slabs,  cylinders  and  spheres)  may  be  used 
to  generate  solutions  for  multidimensional  bodies. 

Use  of  the  product  solution  technique  is  illustrated 
by  Figure  4.2.  The  problem  chosen  to  validate  the  heat 
transfer  portion  of  the  current  model  is  the  "sudden  immer¬ 
sion"  problem  for  a  f inite- length  cylinder  subjected  to 
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Substituting  Equation  4.26  into  Equations  4.22,  4.23  and 
4.24  and  separating  the  variables,  reduces  the  two-dimensional 
problem  into  two  one-dimensional  problems: 
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Thus  Equation  4.26  is  the  exact  solution  to  the  sudden 
immersion  problem  of  a  right  circular  cylinder. 

5 .  Validation  Problem 

The  validation  problem  is  Example  4.9  in  [Ref.  25]. 

The  problem  is  stated  here  for  the  convenience  of  the  reader. 

The  short  cylinder  in  Figure  4.3  is  initially  at  40°  C 

and  then  plunged  into  a  fluid  with  h  =  300  W/m-K  and  T  =  200°  C. 

-9  2 

The  material  is  bronze,  k  =  26  W/m-K  and  a  =  8.6-10  m/s. 
Find  the  temperature  at  the  center  of  the  cylinder  after 
5  minutes. 
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Solution:  176°  C  (348.8°  F).  The  details  of  the  solution 

are  found  on  pages  186-187  of  [Ref.  25].  The  product  solution 

technique  was  utilized  in  obtaining  a  5-term  approximation 

2 

of  the  series  solution  for  Fourier  moduli  (F  =  it/L  , 
equivalent  to  a  nondimensionalized  time)  less  than  0.2.  For 
Fourier  moduli  greater  than  0.2,  a  two-term  approximation  of 
the  exact  solution  was  used.  Heisler  [Ref.  26]  points  out 
that  a  one-term  approximation  of  the  exact  series  solution 
yields  an  accuracy  to  within  1%  of  the  exact  value  for 
Fourier  moduli  greater  than  0.2.  Various  values  (5,  10,  20 
and  30  seconds)  corresponding  to  values  of  Fourier  moduli 
less  than  0.2  were  used  in  comparing  model  solutions  with 
the  exact  (five-term  approximations  to  the  exact)  solution. 
Values  (3,  5,  7  and  9  minutes)  corresponding  to  Fourier 
moduli  greater  than  0.2  were  used  to  compare  model  solutions 
with  the  exact  (one-term  approximations  with  1%  accuracy) 
solutions.  Figures  4.4  through  4.8  show:  (1)  the  effects 
of  time  on  the  solution  of  various  grids,  and  (2)  effects 
of  grid  refinement  at  equal  values  of  time.  Values  plotted 
on  the  ordinate  scales  are  deviation  (percentage  error)  from 
the  exact  solutions  versus  time.  Observation  yields  that  for 
small  values  of  time,  grid  refinement  obtains  more  accurate 
results  and  obtains  a  faster  convergence  to  steady  state. 

6 .  Testing  the  Model's  Ability  to  Accept  and  Reject 
Heat 

An  equivalent  simple  test  in  model  validation  is  the 
model's  ability  to  accept  heat  from  a  boundary  flux  and 
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reject  it  at  a  later  time  to  its  environment  via  convection. 

z  2 

The  applied  heat  flux  at  —  =0  is  1500  Btu/ft  -hr. 

o 

Figures  4.9-4.12  illustrate  system  response.  Figure  4.13 
is  the  input  data  set  used. 

7 .  The  One-Dimensional  Problem 

Another  step  in  the  validation  of  the  two-dimensional 
problem  is  an  examination  of  a  one-dimensional  problem.  One 
might  infer  that  a  two-dimensional  model  includes  the  one- 
aimensional  model  as  a  subset.  This  was  demonstrated  in  a 
separate  test.  A  one-dimensional  problem  may  be  imposed  on 
the  model  in  two  ways.  One  way  is  to  make  the  length  to 
diameter  ratio  very  large.  (No  restrictions  on  excitations 
is  implied. ) 

Another  method  for  eliciting  one-dimensional  behavior 
is  to  "excite  the  system  in  a  one-dimensional  fashion," 
i.e.,  apply  boundary  conditions  at  —  =  0.0  and  —  =  1.0 

2  Z 

o  o 

independent  of  radius  and  insulate  radial  boundaries  at 

=  0.0  and  —  =  1.0.  The  insulated  boundary  at  —  =  0.0 
o  o  ro 

arises  from  an  axisymmetric  formulation  of  the  problem.  The 

zero  gradient  (completely  insulated)  boundary  at  —  =  1.0 

o 

corresponds  to  a  completely  (all  four  variables)  impermeable 
vertical  boundary.  For  this  problem,  an  adiabatic  restriction 
is  sufficient.  In  general,  an  adiabatic  condition  means  an 
insulation  of  heat.  In  this  problem,  the  implication  is  far 
greater  because  of  the  coupling  between  system  fields.  An 
adiabatic  restriction  at  —  =  1.0  prohibits  non-zero  radial 
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pressure  gradients.  Unconstrained  radial  pressure  gradients 

at  —  =  1.0  imply  radial  convection  of  air  enthalpy  and  in  a 
o 

non-isothermal  fluid,  this  convection  undermines  the  initial 

adiabatic  assumption  at  —  =  1.0.  Oxygen  gradients  are  also 

o 

negated  by  an  adiabatic  condition  since  oxygen  convection 
by  air  must  occur  under  the  presence  of  nonzero  pressure 

IT 

gradients  at  —  =  1.0. 

o 

The  two  methods  described  above  distinguish  one-dimen¬ 
sional  behavior  arising  from  geometrical  conditions  and  one¬ 
dimensional  behavior  resulting  from  restrictions  on  excitations. 
Figures  4.14-4.16  depict  model  one-dimensional  solutions  to 
"one  dimensional  excitations."  Figure  4.17  is  the  input  data 
set  used. 
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Figure  1  is  not  available  per  Mr 
David  Salinas,  NPS 


Figure  4.1  Classical  Transient  1-D 
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Infinite  ‘  ] 

cylinder  _ 


Fluid:  h0,  T0 
All  sides: 


Short  cylinder 
formed  by 
intersection  of 
infinite  cylinder 
and  two  planes 


Product  solution: 

0U,  r,  t)-P(x.  t)  C  (r,  f) 


Figure  4.2  Illustration  of  the  Product  Technique 
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PERCENT  ERROR 


Figure  4.5  Grid  Effects  at  Point  B,  (0,1) 


PERCENT  ERROR 


—  O*  «n 


TIMET  (SECONDS) 


Figure  4.6  Grid  Effects  at  Point  C,  (1,1) 
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2 

F  in  the  first  minute)  at  the  —  =  0.0  boundary.  Due  to 

o 

the  steep  gradients,  further  system  behavior  could  not  be 
analyzed  since  output  data  consisted  of  a  single  output 
(T  =  5-20  S.,  problem  time,  depending  on  porosity  value) 
in  15  minutes  of  CPU  time. 

E.  CONCLUSIONS 

The  combustion  analysis  program  (CAP)  is  a  viable  tool 
for  the  analysis  of  heat  transfer  in  porous  media.  The 
model  constructed  provides  the  user  with  the  flexibility  to 
solve  problems  with  or  without  combustion.  The  role  of 
the  permeability  model  must  not  be  underestimated.  For  it 
is  the  physical  assumptions  in  this  aspect  of  the  model  that 
govern  the  flow,  heat  transfer  and  in  the  end  the  evolution 
of  the  combustion  problem. 

F.  RECOMMENDATIONS 

Follow-on  work  is  recommended  in  the  following  areas: 

1.  Develop  a  restart  capability  for  the  program  which 
will  allow  for  a  problem  to  begin  at  any  point  in 
time  with  initial  conditions  and  rate  information 
being  identical  to  previous  timestep  values. 

2.  Nondimensional ize  equations  in  terms  of  other  system 
quantities  in  addition  to  spatial  dimensions. 

3.  Flex  the  model  under  other  system  parameters  to  examine 
effects  on  system  behavior. 
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0.000417  ft.  (0.005  in.).  Values  of  d  (ft.)  equal  to 
0.000417,  0.000411,  0.000396,  0.000385  and  0.000370  were 
used  to  give  values  of  porosity  of  0.476,  0.5,  0.55,  0.60, 
and  0.65,  respectively.  The  results  are  graphically  pre¬ 
sented  in  Figures  5.1--5.4.  Figures  5.1  and  5.2  show  the 
carbon  temperature  profile  and  the  oxygen  concentration 

profiles  for  oorositv  o.  =  0.476  (D  =  0.000417  ft.  =  0.005 
*•  *  -  -  1 

in.).  In  all  cases,  excitations  for  these  runs  are  1500 

2  z  2 

Btu/ft  -hr  (in)  at  —  =  0.0  and  50  Btu/ft  -hr  (out)  at 

o 

IT  Z 

—  =  1.0;  the  oressure  is  14.65  psi  at  —  =  0.0  and  14.55 
r  -  z 

°  z  ° 

at  —  =  1.0;  the  cylinder  length  to  diameter  ratio  is  0.5 

o 

(length  =  1.0  ft.).  Figures  5.1  and  5.2  represent  typical 
graphical  results  for  carbon  temperature  and  oxygen  concen¬ 
tration  profiles.  Figures  5.3  and  5.4  are  a  tabular  summary 
of  the  results  of  the  carbon  temperature  profile  and  the 
oxygen  concentration  profiles,  respectively,  for  varying 
porosities.  Profiles  with  higher  values  of  initial  porosity 
are  observed  to  have  accelerated  development  of  temperature 
profiles  (i.e.,  as  the  fuel  diameters  decrease  (increasing 
porosity) ,  the  carbon  temperature  responds  at  a  faster  rate 
to  system  excitations) .  The  oxygen  concentration  response 
is  as  expected  (i.e.  as  reaction  rate  goes  up,  the  oxygen 
concentration  goes  down).  In  the  runs  made  for  porosity 
values  of  0.70,  0.75,  0.85  and  0.90,  numerical  difficulties 
were  encountered.  The  temperature  profiles  were  observed  to 
increase  sharply  compared  to  lower  values  of  porosity  (300-400° 
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problems  discussed  here  the  transport  crooerties  ,  h  , 
m/u ,  Ce  and  the  Arrhenius  reaction  rate  term  are  :cvern;rm 
factors  in  the  respective  field  equations  in  which  they 
appear.  As  previously  discussed,  the  net  balance  between 
heat  transfer  rates  and  heat  generation  dictate  whether  a 
reaction  will  proceed  to  extinguishment  or  combustion. 

There  are  two  time  constants  to  be  considered  in  a  combustion 
problem  of  this  type:  a  time  constant  for  the  Arrhenius 
reaction  and  a  time  constant  associated  with  momentum 
transport.  The  momentum  time  constant  affects  the  rate  at 
which  heat  is  convected  (removed)  out  of  a  differential 
volume.  The  reaction  time  constant  affects  the  rate  at 
which  heat  is  generated  (added  into)  within  a  differential 
vo 1 ume . 

3.  POROSITY  ANALYSIS 

In  this  problem  there  are  many  parameters  one  might  vary 
in  order  to  examine  resultant  system  behavior.  Due  to  time 
limitations  this  investigation  examines  the  effect  of  one 
parameter,  porosity,  on  system  behavior.  The  observed 
effect  of  porosity  values  ranging  from  0.476  to  0.90  is 
dis  '  ^d.  Of  the  nine  runs  attempted,  only  five  had  suffi¬ 
cient  ut  data  that  allowed  for  comparative  analyses. 

Employi  'uation  3.1  for  the  porosity  associated  with 
spherical  particles,  the  carbon  diameter,  d,  was  varied  to 
achieve  various  values  of  porosity  while  holding  the  parti¬ 
cle  center- to-center  distance  D  (Figure  3.1),  fixed  at 
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Similarly,  in  place  of  Equation  3.44,  a  zero  flux  boundary 

condition  at  =  1.0  could  have  been  imposed  on  the  oxygen 
o 

concentration.  The  zero  flux  condition  implies  an  impermeable 

boundary  with  respect  to  oxygen  concentration  fluxes.  This 

boundary  condition  would  apply  to  a  very  long  cylinder  (i.e., 

a  regularity  condition).  In  this  investigation,  a  Cauchy 

(convective  flux)  boundary  condition  on  the  oxygen  concentra- 

z 

tion  is  imposed  at  —  =  0.0  and  1.0.  As  oxygen  is  locally 

o 

consumed  in  the  interior  of  the  system,  oxygen  gradients  are 
created.  According  to  Fick's  Law  of  Diffusion,  a  depleted 
oxygen  region  is  replenished  by  the  diffusion  of  oxygen  from 
regions  of  relatively  high  concentration  to  regions  of  low 
concentration.  Thus  a  convective  flux  boundary  condition 
on  oxygen  causes  the  depletion  of  oxygen  concentration  to  be 
retarded  through  the  diffusion  mechanism. 

B.  EXCITATIONS 

The  effects  of  high  and  low  values  of  the  heat  flux 
2 

applied  at  —  =  0.0,  are  predictable.  High  values  of  heat 
o 

flux  lead  to  an  accelerated  development  of  both  the  carbon 

and  air  temperature  profiles  and  0^  concentration  depletion 

within  the  system.  For  high  values  of  heat  flux  at  —  =  0.0, 

zo 

the  numerical  integration  scheme  eventually  slows  noticeably 
as  a  result  of  steep  gradients  observed  at  this  boundary. 

C.  PHYSICAL  CHARACTERISTICS  OF  THE  SYSTEM 

The  physical  characteristics  of  the  system  generate  a 
significant  effect  in  field  problem  solutions.  In  the 
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V.  DISCUSSION  AND  CONCLUSIONS 

Field  problem  behavior  depends  on  three  major  factors r 
(1)  Boundary  conditions,  (2)  System  excitation,  and  (3)  The 
physical  characteristics  of  the  system. 

Time  limitations  did  not  permit  this  investigator  to 
conduct  an  exhaustive  analysis  of  each  of  these  factors. 

Some  preliminary  analyses  were  performed  to  obtain  some 
understanding  of  the  effect  of  boundary  conditions,  and 
excitation  on  system  behavior. 

There  are  a  variety  of  physical  parameters  which  govern 
system  response,  such  as  permeability,  porosity,  and  the 
cylinder  length-to-diameter  ratio.  Here  only  a  brief 
investigation  of  the  effect  of  porosity  on  system  behavior 
was  undertaken,  and  is  reported. 

A.  BOUNDARY  CONDITIONS 

The  boundary  conditions  used  in  the  present  investiga¬ 
tion  were  presented  in  Chapter  III.  Other  boundary  conditions 

are  possible.  For  example,  if  in  place  of  Eqs .  3.32  and  3.36, 

2 

insulated  boundary  condition  at  —  =  1.0  is  imposed  on  the 

o 

carbon  and  air  temperatures,  then  preliminary  results  indi¬ 
cate  temperature  response  is  higher  for  equal  values  of  time. 
This  behavior  is  due  to  the  buildup  and  storage  of  energy 
within  the  system  associated  with  an  insulated  boundary. 


an 
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Figure  4.15  One  Dimensional  Problem  (Carbon 
Temperature  Profile) 
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One  Dimensional  Problem  (Carbon 
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Figure  4.11  Model  Accepting  Heat  (Oxygen 
Concentration  Profile) 
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4.  Consider  an  internal  ( particle-to-particle)  radiative 
heat  transfer  analysis . 

5.  Consider  other  fuels  and  more  detailed  chemical 
kinetics  chains. 
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Figure  5.3  Maximum  Carbon  Temperature  Summary  Porosity: 
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Figure  5.4  Oxygen  Concentration  Profile  Summary 

Porosity:  p,  =  .476,  p2  =  0.5,  p,  =  0.55, 
p^  =  0.60,  p^  =  0.65.  Time:  t^  =  0.1 
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APPENDIX  A 


FORMULATION  OF  FIELD  EQUATIONS 

A.  PRESSURE  DISTRIBUTION  EQUATION 

Darcy's  law  for  two-dimensional  flow  is, 


Q 


E(3P  r  + 

U  or  o  z 


p  pa  g)  z) 

a. 


(A.  1) 


where  Qr  and  Qz  are  expressed  as  follows, 


Q 


r 


m  _3p 
la  3r 


m  3P  . 

-  (t—  ~  pp  q) 
y  o  z  a 


(A. 2) 


(A. 3) 


Invoking  the  Dupuit-Forcheimer  relation,  and  solving  for  the 
pore  velocity  components  u  (radial  velocity)  and  v  (axial 
velocity),  Equation  A.l  becomes. 


-m  3P 
pu  3r 


(A. 4) 


v 


—  m  ,  3  P 
pu  9z 


P  Pag) 


(A. 5) 


The  continuity  relation  (derived  in  Appendix  B)  is, 


D (ppa) 
Dt 


+  P  Pa  Div  v 


0 


(A. 6) 
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Substituting  Equation  A.  2  and  Equation  A. 3  into  Equation  A. 6 
yields , 


3p 


a  g  ,  m  ,  3p  ~  ,  ,  3P  .  ~ . 

—  =  p  pa?-  W(^  r  '  P  °ag)z) 


+  ;  +  -  PV>;)-?poa 


(A. 7) 


Expanding  terms,  Equation  A. 7  becomes, 


32P  +  32P  +  ,_1_  _^a  +  1  3m 
3?  7?  pa  3r  m  3r 


1  3U  .  1, 

U  3r  r  3r 


•  /_L  3pa  +  I 

p  3z  m  3z 

a 


1  ,3P  . 

U  32>  (1T  +  P  pag): 


3Po 

u  —  ( A.  8  ) 


pm  3 1 
3 


Equation  A. 7  with  associated  boundary  conditions  (presented 
in  Section  II. B)  is  cast  into  a  finite  element  formulation 
and  becomes  one  of  four  field  equations.  Pressure  gradient 
information  is  then  substituted  into  Equation  A. 4  and  Equation 
A. 5  to  obtain  the  pore  velocity  components. 


B.  POROUS  SOLID  HEAT  TRANSFER  EQUATION 

In  performing  energy  balances  on  both  the  porous  solid 
and  on  the  air,  a  differential  volume  of  porous  medium  may 
be  partitioned  into  respective  volumes,  dVg  =  (l-p)dV  for 
the  solid,  and  dV  =  pdv  for  the  air  (shown  in  Figure  A.l) 
The  convention  used  for  the  energy  balance  of  an  arbitrary 
differential  volume,  dV,  is. 
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Heat  into  DV  +  Heat  generation 


Heat  out  of  DV 


+  Increase  in  internal  energy  (A. 9) 

The  heat  transfer  mechanisms  considered  for  the  carbon 
are  conduction,  radiation  heat  transfer  between  particles, 
convection  heat  transfer  from  the  particles  to  the  air, 
and  heat  generation.  Applying  the  above  convention,  the 
energy  balance  on  a  differential  volume  of  porous  solid  is 
(invoking  Tayler  Series  expansions  and  neglecting  higher 
order  terms) , 

1-  7  1F((1-P>r<Jcond,r>  -  T7((1-P>qco„d,z> 

-  IfklU-plrq^p  -  |j<(l-p)qradi2)]dV 

-  qcor,vdA'  +  qgendA’  =  (1-P)qintdV  <A'10 

In  vector  form,  Equation  A. 10  becomes. 


- 7 ■ ( q  ,  +  q  , ) dV  -  q  dA '  +  q  dA ' 
^cond  ^rad  conv  ^gen 


(1-p,qintdV 


(a.  n; 


Substituting  the  following  expressions  into  Equation  A. 10, 


* 


rad 

=  -  k  7T 

r  c 

Fourier  analogy 
( radiative) 

(A. 13) 

conv 

=  h  (T  -  T  ) 
c  a 

Newton's  Cooling  Law 

(A. 14) 

gen 

=  R 

g 

Heat  generation 

(A. 15) 

int 

=  .  c 

Kc  c  at 

Internal  energy 

(A. 16) 

yields , 


3  T 

ii|_(ra-p)0<e)  3r,  .  8z 


c)  +  f-( (1-P) (k 


3T 

)  — 

e  3  z 


•)  )dV 


h ( T  -  T  ) dA '  +  R  dA ' 
c  a  g 


oT 


-  (1'p)pc  Ccft 


dv 


(A. 17) 


Dividing  through  by  dV,  and  defining  dA'/dV  as  Z,  the  specific 
internal  area  (i.e.,  surface  area  per  unit  volume),  Equation 
A. 17  becomes, 


7*  (  ( 1-p)  (k  ) VT  )  -  hZ (T  -  T  )  +  RZ  =  (l-p)p  C  £  (A. 18) 

SC  Cd.  y  C  C  d  l. 

The  expressions  used  to  obtain  values  of  the  properties  and 
parameters  in  Equation  A. 18  are  presented  in  Section  III.E. 


C.  AIR  HEAT  TRANSFER  EQUATION 

The  formulation  of  the  air  heat  transfer  equation  begins 
with  the  general  two-dimensional  energy  balance  equation, 


★ 

The  difficulty  in  obtaining  an  expression  for  kr  is 
addressed  in  Section  III.E.  Throughout  this  work,  one  may 
keep  in  mind  that  ke  should  really  be  (ke  +kr)  to  account 
for  conduction  as  well  as  radiation  within  the  porous  medium. 
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(A. 19) 


p  p  D  (U  +  K) 

cl 


=  -7 • pq  -  V-pPV  -  (7*(px*V)) 


+  hZ(T  -T  )  +  p  pa  (V-g) 
c  d  a. 


where  U  is  the  internal  energy,  K  is  kinetic  energy,  t  is 

the  dissipation  function,  and  g  is  the  gravitational  acceleration. 

Expansion  of  Equation  A. 19  yields, 


D  .  1,2  2.  . 

P  °a  Dt  e  +  2(u  +v  )} 


-7-  ( -pk  7T  )  -  7*  (p  P  V) 

a.  a. 


-II  g -  p  t.  .  V  .  +  hZ  (T  -  T  ) 

i  i  3xi  13  3 


+  P  P, v  g 

a 


(A. 20) 


where 


(PT-V)  =  y  y  i p  T.  .V.  =  7  ( px  v  +pT  v 

^  £  t-  9x^  r  i]  ]  ^  rr  r  *rz  z 

+  px  V  +  px  V  ) 
^  zr  r  c  zz  z 


|j(PTrrU  +PTrzV)  +  fj(PTzrU  +PTZZV)  (A. 2i: 


And  so  Equation  A. 19  becomes, 


rDe  0,1,2  2 

Dt  +  dt  2  u  +v  1 


*  V  (pkaVTa)  -  ^(pirrU  +P  TrIv) 


“  Iz(pTzru  + pTzzv)  ‘  (VpP‘ v  + Ppv’ v) 


+  hZ(T  -T  )  +  p  p  vg  •  (A. 22) 

c  a  a 
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As  in  the  porous  solid  heat  transfer  equation,  the  time  and 
position  dependent  porosity  appears  inside  the  differential. 
Expanding  terms,  the  air  energy  equation  becomes, 


De 


PP. 


a  D  ,  2  2, 

p  +  ~T~  Dt  u  +v 


V*  (p  k  9T  )  -  [7  p  P-V  +  p  P  7-V] 
S3 


-  ( pr  u  +  pi  v) 

or  r  rr  r  rz 


-  ■£—  (pi  u  +  pi  v)  +  hZ  (t  -T  ) 
o z  ^  rz  r  zz  c  a 


+  P  Pav  9 


Consider  the  momentum  equation  for  the  r-direction, 


R  Momentum 


ft tp  cau] 


1,3,  2 .  3  .  . 

r  ("  3?(r  ppa  u  >  '  3?(p  pa  UV)  3r 


3pP 


-  (—  I—  (p  r  x  )  - 
r  3r  r  rr 


pi 


89  ,  3 


+  ^-(p  x  )  ) 
r  3z  r  rz 


Consider  the  momentum  equation  for  the  z-direction, 


Z  Momentum 


It(p3aV) 


1,  3  ,  . 

—  -  r—  r  p  p  u  v 
r  or  ^  a 


tj(p  Cav‘)  -  fjlpp) 


•  '?  f?<p  rrz  1  +  I V(p  Tzzn  +  p  pa9 
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(A. 23) 


(A. 24) 


(A. 25) 


and  the  continuity  equation, 


Equation  A. 27  becomes 


Du 

P  °aU  Dt  = 


3  pP  ,13,  , 

-u  ^ - u  (—  T— (p  rr 

a  r  r  or  rr 


PT  ,  - 


T~  +  T7(P  ~  rz  ]  ] 


(A. 30) 


Multiplying  Equation  A. 28  by  v  and  noting  that, 


Dv 

D  0  V  — 

^  a  Dt 


j  v  3  v  2  3  v 

pov^rr+pouv  —  +  p-v  — - 
r  a  ot  a  or  a  oz 


{A. 31) 


Equation  A. 28  becomes 


Dv 

p  “aV  Dt  = 


-  v 


)pP 


3  Z 


V(F  fc'P  rTn"  +  fl(pTzz)  1  +  P  "a''  g 


(A. 32) 


Expanding  Equation  A. 30  yields, 


p  °a  U 


Du 

Dt 


u 


3pP 

3r 


3pi 


rr 


3r 


u 

-—pi 

r  c  rr 


upxai 


-  u 


4-(px  : 

oz  c  rz 


(A. 33) 


Expanding  Equation  A. 32  yields. 


Dv 

P  °aV  Dt 


3pP  3(pxrz: 

‘V  3T  “  V  3r 


v 

7:  P 


rz  -  v  H'P'zz1  +  p  •'aV  9 


(A. 34) 


The  energy  Equation  A. 23,  after  substituting  and  expanding 

terms,  is, 
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The  viscous  dissipation  terms  in  Expression  A. 36  are 
neglected  because  the  fluid  is  a  gas  flowing  at  low  velocity 
(see  [Ref.  9]).  Therefore,  the  energy  equation  for  the  air 
in  the  porous  medium  is, 
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AD-A151  965  TfclO-D  HEAT  TRANSFER  THROUGH  POROUS  MEDIA  WITH  HEAT  2/2 

GENERATIONS)  NAVAL  POSTGRADUATE  SCHOOL  MONTEREY  CA 
B  MARTINEZ  JUN  84 


UNCLASSIFIED 


F/G  28/2 


NL 


MICROCOPY  RESOLUTION  TEST  CHART 


(A. 37) 


P  “a  Dt 


De  _ 


;pk  7T  )  -  pPDiv  V  +  hZ(T  -T  ) 

3  3  C  3 


With  specific  enthalpy  for  a  gas  defined  by, 


ii  =  e  +  — 


(A. 38) 


De/Dt  can  be  expressed  as, 


De 

Dt 


DH 

Dt 


1  DP 


Do 


Dt 


2 
'  a 


Dt 


[A. 39) 


where  k  is  the  specific  enthalpy,  and  e  is  the  specific 


energy  (internal  plus  kinetic).  Multiplying  the  continuity 


equation  through  by  P/(pp  )  obtains, 

3. 


-£l!t(ppa)  +-Z7?fc(pparu)  +  Jrfe(ppav) 


=  0  .  (A. 40) 


PP. 


PP. 


PP. 


Expanding  yields. 


PP 


JL  +  o  l£) 
2XfcJ  It  Pa  dt’ 


■(P 


P  3  ,  , 

— 7u37tE,pa’ 
?pa 


PP 


a  3 


PP. 


2  r  3r 


(ru) 


P  v  3  ,  , 

-7  3pe°a; 
PP. 


P 


3  v 


p  pa37 


(A. 41) 


In  vector  notation,  Equation  A. 41  becomes, 


(A. 42) 
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'p  T 


-  dP 

P 


(A. 48) 


and  recalling  that, 


_  1 

p  3T  p 


(A. 49) 


the  equation  of  state  for  a  perfect  gas, 


P  =  p  RT 


(A. 50) 


simplifies  Equation  A. 49  such  that, 


3  =  - 

RT 


(A. 51) 


Thus,  Equation  A.  4  8  becomes, 


ds  =  C  ~  -  dP 
p  T  pT 


(A. 52) 


Substituting  Equation  A. 52  into  Equation  A. 47  and  cancelling 
terms  yields. 


C  dT 
P 


(A. 53! 


a  Dt 


(A. 54) 


98 


Substituting  Equation  A. 54  into  Equation  A. 45  yields. 


9T 


P  C 


a  a  3t 


V- ( pk  VT  )  -  pp  C  (V- V) T 

d  cl  d  d  d 


+  hZ  (T„  -  TJ  +  ° 


c  a 


(A. 55) 


Making  the  assumption  that  pP  changes  very  little  with  time 
[Ref.  9],  (this  assumption  was  subsequently  confirmed  by  the 
model) ,  i . e . , 


DpP 

Dt 


(V-V)pP  =  u 


or 


+  v 


a  (pp) 

9z 


(A. 56) 


The  final  air  heat  transfer  (energy  conservation)  equation 
is, 


p  °aCa 


7-  (pkaVTa)  -  p  OaCa(V-V)Ta  +  u  +  v 

+  hZ (T  -  T  ) 
c  s 


iE£ 

dz 


(A. 57) 


In  vector  form.  Equation  A. 57  becomes, 

ppaCa  ~  5-(pkaVTa)  -  ppaCa(5-V)Ta 

+  (V'V)pP  +  hZ(T  -T  )  (A. 58) 


The  expressions  used  to  obtain  the  properties  and  parameters 
in  the  coefficients  of  Equation  A. 57  are  presented  in 
Section  III.E. 
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D.  OXYGEN  MOLECULE  DIFFUSION  EQUATION 

The  final  consideration  in  formulating  the  field  equations 
for  the  model  is  the  transport  of  oxygen  molecules.  The 
oxygen  molecule  transport  equation  is  obtained  by  a  conser¬ 
vation  of  species  balance  on  the  differential  volume  of  air, 
dV  =  pdV.  The  convention  used  for  the  species  balance  into 
a  differential  volume  is, 


C>2  in  =  C>2  out  +  C>2  consumed  +  O2  accumulated  (A.  59) 

The  transport  mechanisms  considered  were  diffusion  due  to 
concentration  gradients  (Fick's  law),  convection,  and  air 
consumption  by  combustion.  The  species  balance  on  oxygen 
becomes , 


pm , . - -dA  +  pm,  .--dA  +  pm  dA  ,  +  pm  dA  j 

^  diff  ^  diff  r  conv  r  conv  ' 

r  z  1  r  z 


pm ,  . --dA  +  pm,  .--dA  +  pm  dA, 

^  diff  ,  .  r  diff  r  conv  , 

r+dr  'z+dz  r  +  dr 


+  pm  dA 
r  conv 


z+dz 


+  m  dA 1  +  pm  dV 
cons  acc 


[A. 60) 


Representing  terms  on  the  right  side  by  Taylor  series  expan¬ 
sions  (neglecting  higher  order  terms),  i.e.. 


pm^dA 


ei+dei 


=  pm,  dA 


3ei 


(pm 


kdA)d£i 


(A. 61) 
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where , 


k  =  conv,  =  r  or  z  dAf  =  r  d0  dZ 


and  dAz  =  r  dr  d0 


Then  Equation  A. 61  becomes, 

3  •  3 

t— (p  mdA  )dr  =  (r  pu  i  d0  dz)  dr 

JI  L  JL 

=  -p  |^r(r  pu  4>)  r  de  dZ  dr 

and 

3  •  3 

(p  mdAz)dz  =  v$  r  dr  dQ)  dz 

£ 

=  -^-(p  v  i)  r  d9  dz  dr 


(A. 62) 

(A. 63) 

(A. 64) 


Substituting  Equations  A. 64  and  A. 65  into  Equation  A. 60, 
cancelling  terms  and  rearranging.  Equation  A. 61  becomes, 


-3?(P(mdiff 


conv 


dr  -  3I(p(mdiff 


m  )  dA  ; 
conv  z 


m  dA ' 
cons 


=  pm  dV 
acc 


(A. 66) 


Substituting  the  following  expressions  into  Equation  A. 66, 
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POROUS  MEDIUM 


AIR 


SOLID 


Figure  A.l  Separating  A  Differential  Volume  of 

Porous  Medium  into  Respective  Volumes 
of  Solid  and  Air 


(A. 67) 


m 


diff 


Ve  7  0 


m 

conv,r 


ug , 


m 

conv,  z 


vo 


m  =  R_ 

cons  C>2 

•  _  3(f) 

m  — 

acc  3t 


(A. 68) 


(A. 69) 


(A. 70) 


dividing  both  sides  by  dV,  and  letting  dA'/dV  equal  the 
specific  internal  area,  z,  the  oxygen  molecule  diffusion 
equation  becomes, 


7 *  (pPe7$) 


v*  (p  o  v) 


p  H 

3t 


(A. 71) 


The  methods  and  expressions  for  obtaining  the  properties  and 
parameters  in  the  coefficients  of  Equation  A. 71  are  presented 


in  Section  m.G. 


APPENDIX  B 


AUXILIARY  EQUATION  FORMULATION 

A.  CONTINUITY  EQUATION 

The  continuity  equation  for  a  fluid  in  a  porous  medium 
is  expressed  by, 

D(po  ) 

- — -  +  p  f>a  Div  V  =  0  (B.l) 

or  in  equivalent  form, 

hripoj  +  (V-V)pc  +  pp  (7  •  V)  =  0  (B .  2 ) 

j  u  d  ct  d 


and 


jpp 

+  u  — =; -  +  V 

or 


0(P "a} 

oZ 


+  Ppa(? 


9  (ru) 
9r 


+ 


( B  .  3  ) 


B.  LAGRANGE  POLYNOMIAL  APPROXIMATIONS  FOR  THERMAL  PROPERTIES 
Relations  for  calculating  the  dynamic  viscosity,  thermal 
conductivity,  and  specific  heat  at  constant  pressure  of  air 
at  different  temperatures  were  required.  Second  order 
Lagrange  polynomial  fits  to  empirical  data  provide  a  simple 
method  fcr  obtaining  the  relations  required.  Vatikiotis 
[Ref.  9]  gives  the  details  for  arriving  at  the  resulting  set 
of  polynomials, 


104 


■< — w 


V  ■  = 


*  * .  .  ~  ’  v r* r ^  '  t  r ». 1 J9  f  y*  j  «  r» ■  ■ ', » ?  >■  ■  ■:»  y  ■  ■■  »■ 


=  -3.308  xio'9  T2  +  4.633  x  10~5  T  +  4.427  x 10  2  (B.12) 

Pi  cl 


k.  =  -2.608  xio'10  T2  +  1.930  xio'5  T  +  1.361  x  10  2  (B.13) 

cl  ^  <3. 


Ca  =  -1.293  x 10"9  T2  +  2.758  x 10  5  +  .238  (B.14) 


Each  expression  obtains  property  values  within  two  percent 
of  the  data  presented  in  [Ref.  9]  for  temperatures  up  to 
3000  degrees  Fahrenheit. 


C.  CARBON  PARTICLE  SURFACE  RECESSION 

The  following  analysis  of  particle  diameter  consumption 
assumes  that  the  fuel  particle  surface  recedes  uniformly. 

This  assumption  is  reasonable  if  the  particle  diameter  is 
small  in  comparison  to  the  cylinder  geometry,  i.e.,  negligi¬ 
ble  boundary  effects.  In  addition,  since  the  velocities  are 
low,  the  hydrodynamic  effects  on  the  uniformity  of  the 
particle  surface  recession  are  negligible.  The  analysis  also 
assumes  there  are  no  significant  thermal  gradients  within 
the  particle,  i.e.,  an  isothermal  carbon  particle.  This  is 
also  reasonable  for  the  small  particles  examined  in  this 
analysis.  A  mass  balance  must  be  performed  on  the  particle 
and  equated  with  the  reaction  rate.  This  equivalence  yields 
the  following  expression, 


dm 

dt 


R  Z 
c 


(B. 15) 
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APPENDIX  C 


GALERKIN  FEM  FORMULATION 

A.  FINITE  ELEMENT  METHOD 

The  solution  of  the  system  of  coupled,  nonlinear  partial 
differential  equations  given  by  Equations  3.9,  3.22,  3.26 
and  3.27  —  subject  to  boundary  and  initial  conditions,  was 
obtained  by  a  Galerkin  formulation  of  the  Finite  Element 
Method . 

1 .  Galerkin  Formulation 

A  Galerkin  formulation  of  the  Finite  Element  Method 

was  used  to  obtain  solutions  of  the  porous  solid  and  air 

energy  equations,  the  oxygen  diffusion  equation,  and  the 

continuity  (pressure--Darcy ' s  law)  equation.  A  convenient 

form  of  Equations  3.9,  3.22,  3.26  and  3.27  was  used  in  the 

formulation  where  the  spacial  coordinates,  r  and  z  were 

nondimensional ized  by  e  =  z/z  and  n  =  r/r  . 

o  o 

The  closed  domain  defined  in  (r,z)  space  by  (0,0), 

(0,1),  (1,1),  and  (1,0)  was  partitioned  into  NEL 

(2*  (NRNP-1)  MNZNP-1)  )  contiguous  area  elements  obtained 

from  a  NSNP  model.  NSNP,  NRNP  and  NZNP  are  the  number  of 

system  nodal  points,  radial  points,  and  axial  points, 

respectively.  The  four  field  variables  T  ,  T  ,  P  and 

c  a 

were  approximated  by. 


107 


121 


term  is  rearranged  as  follows, 


AR"  n  lexp{-E/RuTc)  , 


‘V1 


(C. 44) 


Letting 


:c.  45; 


_T . 

z  ;4 


[C. 46) 


and  invoking  natural  coordinates  [Ref.  28],  an  elementally 
averaged  contribution  results  in  the  following  area  integral, 


c  //  5 

Z  A 

e 


44  dA 


(C. 47; 


Expanding  yields, 


-  ,  ,  T 

c  =  c  Ij 


•l^j  !  ’2 ^ j  :  ’ 3  ’ j 


1  i  4  ] 


:c. 48) 


1  =  1,3,  j  =  1,3 


A  final  explicit  expansion  of  the  integral  looks  like, 
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The  reaction  rate  term  may  be  treated  in  various  ways, 

a.  Treatment  as  an  Excitation  (Force  Term) 

This  treatment  is  the  easiest  of  all  schemes. 

The  term  is  merely  evaluated  at  the  last  time  step  and  is 
assumed  to  be  constant  over  the  next  integration  time  step, 
i  .  e .  , 


R  =  {A  rJJ  exp  (E/R  T  )  }  (C.42) 

The  superscript  *  indicates  evaluation  occurring  at  the 
previous  time  step.  The  value  of  the  term  is  evaluated 
at  the  ith  nodal  point  and  inserted  as  a  factor  in  the 
corresponding  carbon  energy  or  the  oxygen  diffusion  equation 
of  the  ith  nodal  point. 

b.  Linear  Operator  Treatment  of  C>2  Concentration 

In  order  to  realize  an  improvement  over  the  first 
treatment,  one  may  retain  a  portion  of  the  reaction  expression 
as  an  operator  by  making  the  following  rearrangement: 

R  =  (A  R™  <J>n_1exp(-E/R  T  )  }*•  [<£]  (C.43) 

g  u  c 

where  [f]  denotes  a  spatial  operator  treatment  of  the  response 
variable . 

c.  Temperature  and  0_  Concentration  Bilinear 
Operator  Treatment  of  the  Reaction  Rate  Term 

The  bilinear  operator  treatment  of  the  reaction 

rate  is  the  present  method  of  treatment.  The  reaction  rate 
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Cauchy  boundary  conditions  is  as  follows. 


<p  k7;<d  E, 


j>  k  ( ;<  -  X, 

ie 


k  — 

Kie  3 


—  — V 

F  “1 

2  1 

*1 

i  i 

X 

2 

L» 

^  0 

scattered  into  3  and  -  K ,  -r-  added  to 

G  2 


F ( 1 )  and  F(2)  at  the  element  nodal  points 


(C. 38) 


For  a  complete  development  of  the  theory  for  incorporating 
boundary  conditions,  see  [Ref.  28]. 

3 .  Treatment  of  the  Reaction  Rate  Term 

An  exponential  reaction  rate  term  appears  in  both 
the  porous  solid  and  oxygen  diffusion  equation.  The  reaction 
expression  is, 


Rc  =  A  exp(-E/RuTc) 


(C. 39) 


In  the  carbon  equation,  Expression  C.39  appears  as. 


r  Z  =  (R  AH  ) Z 

g  c  K 


(C. 40) 


In  the  oxygen  concentration  equation,  Expression  C.39 
appears  as, 


V  = 


Rc£R1)Z 


( C  .  4  1 ) 
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2 .  Implementation  of  Boundary  Conditions 


Having  formulated  the  system  matrices  for  the  field 
equations,  treatment  of  the  boundary  conditions  is  now  dis¬ 
cussed.  Each  field  equation  is  treated  individually. 

a.  Porous  Solid  Energy  (Heat  Transfer)  Equation 

There  are  three  types  of  boundary  conditions  that 
can  occur:  Dirichlet,  Neumann  and  Cauchy  (mixed)  boundary 
conditions.  The  treatment  of  each  is  as  follows.  For 
Dirichlet  boundary  conditions,  one  may  specify  an  equation 
to  be  a  linear  equation  (i.e.,  independent  of  time)  by  placing 
a  1  on  the  diagonal  of  the  system  stiffness  matrix  correspond¬ 
ing  to  the  particular  degree  of  freedom  at  hand  and  setting 
each  time  derivative  matrix  coefficient  for  the  same  equation 
equal  to  zero.  An  alternative  scheme  involves  setting  the 
time  derivative  diagonal  term  equal  to  one  and  setting  all 
stiffness  coefficients  (in  the  same  equation)  equal  to  zero. 

In  this  manner  one  specifies  the  Dirichlet  value  as  an  initial 
condition  whose  residual  is  identically  zero  and  thus  is 
invariant  with  time.  Treatment  of  Neumann  boundary  conditions 
is  as  follows, 

,  ,  ^  Ze 

4  k7xd£;  =  <j>  p  d£  -*■  “  added  to  local  F(l) 

le  le 

and  local  F(2)  at  the  element 
nodal  points  (C.37) 

In  the  above  expression,  ;<  is  a  response  variable  and  g  is 
an  average  value  of  flux  across  an  element.  Treatment  of 
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; j  NJ> )  dA 

A  ^  r 
e 


— [b  .  b  .  ]  5  _ 
4 A  1  3  -elt 


(C. 31) 


If  N  (b)  zdA  = 

A  " Z  Z 
e 


- [c  .  c  .  ]  6  ,  , 

4A  1  ^  -elt 


(C. 32) 


//  NvdA  =  T|[k]0elt,  k  =  1,  i  ji  j 
A 

e  ... 

k  =  2,  1=3 


(C. 33) 


<p  kN  ( b )  d-',  =  $  kN(.)  dr  -  <J>  kN(v)  dz  (C.34) 

■  .  n  ^  -r  r ,  ~  z 

<.  e  ..  e  l  e 

where  a(i,j)  coefficients  of  the  element  matrices  A(3  *3)  are 
given  by, 

a(i,j)  =  [gl  ;  g  =  g(i,j)  (C.35) 


and  the  vector  0e^t  is  given  by, 


A 

-elt 


(C. 36) 


The  derivations  of  these  operators  are  presented  in  Section 
A. 5  of  this  appendix. 
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To  formulate  these  operators,  the  global  linear 
shape  function,  ,  are  defined  on  the  local  level  by. 


(C. 25; 


where  the  natural  coordinates  ^  anc*  C3  are  defined  by 

Figure  C.2, 


s  ri,  i  =  1,2,3 


(C. 26; 


and  A  is  the  area  of  the  elemnet  and  the  A.  are  the  areas 
e  1 

(in  Figure  C.2)  subtended  by  lines  from  a  point  P(r,z) 
inside  the  triangle  to  the  triangle's  vertices  (rj,Zj, 
j  =  1,2,3).  The  local  shape  functions  (i.e.,  the  elements) 
have  the  following  properties, 


ip.  .  =  9  •  • 

1] 


3  "] 


,  i  -  1,4,  j— 1,3 


;c.27) 


S  i(NPj)  =  6ij  ;  i  =  1,3,  j  = 


( C .  2  8 ) 


Having  defined  the  local  shape  functions,  the  ele¬ 
mental  matrix  operations  C.19  through  C.24  are, 


/ /  N(v) r  dA  = 


f[bj39elt 


(C.2 9) 


II  N(l)  dA  = 


6  [cj 1 ?elt 


(C. 30) 
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Implementation  of  the  boundary  conditions  is  presented  in 
Section  2  of  this  appendix.  The  coefficients  in  Equations 
C.15,  C.16,  C.17  and  C.18  are  temperature  dependent  properties 
and  were  taken  as  the  average  values  of  the  properties  over 
an  element.  In  the  limit,  as  the  elements  get  smaller  (i.e., 
NSNP  -  ») ,  the  average  values  of  the  coefficients  converge 
to  the  exact  values  Inspection  of  expressions  C.15,  C.16, 
C.17  and  C.18  yields  the  six  operators, 


JJ  NO)  dA 


(C.19) 


// 

A 


N  ( 'i) )  dA 


(C. 20) 


II  N  O)  dA 

A  r 


( C . 21 ) 


II  N  CP)  dA 
A  ~z  z 


(C . 22) 


II  N:j,dA 

A  ~ 
e 


and 


1 1  NipdA 
A  ~ 


e 


(C .  2 3 ) 


<£  N(l)  d  l 


(C. 24) 


//NRdA--*  ^N(P)  d2  +  //  -2-  N  (P)  dA 

V  p  ;e  Ae 


ff  pJ  9po 

/)  _  N2,p,zdA  *  ^  ^  N  ( P )  r d A 

e 


-  // 


m  _  a 
py  9  z 


N(P) zdA  +  If  -B-  N  PdA 
A  a 


(C. 17) 


I  JNR.dA  = 

A  ~  4> 

e 


^.epPe?(4>)ndi  +  //  pDe*rU)rdA 

A 

e 


pP 

II  _£  N(4>)  dA  -  £  pP  N(M  d£  +  //  pP  N  ( p )  dA 
A  r  ~  r  £e  e'*  n  A  e~z  z 

e  e 


( rpu) 


+  //  puN  (<}) )  dA  +  //  pvN  (4>)  dA  +  //  - — — —  N4>dA 

a  -  IT  .--Z  m  i- 


+  //  (pv)  N<()dA  +  //  R  ZNdA  +  //  pN0dA 

Ao  Z~  Ao  U2  ~  A  ~ 

e  e  e 


(C. 18) 


The  line  integral  terms  in  each  expression  above  are  boundary 


terms  which  permit  incorporat ion  of  natural  boundary  conditions, 
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*fePka?(Ta’ndi  +  [>  pka?r(Ta'rdA 

e 


//  ^  N(T  )  dA  +  JJ  pkaN2(Ta)zdA 
Ae  Ae 


j>  PkaN(Ta)nd£  -  JJ  hZN(Tc-Ta)dA 

A  _ 


te 


1 1  upN (P)  dA  -  // 
Ae  ~  Ae 


u  NPdA 

3r  ~ 


JJ  vpN(P)  dA  -  //  v  NPdA 
A_  ~  Aq 


dA 

z 


//  P»aCa?Ta  ^ 


Adopting  the  convention, 


NSNP 


=  N  9  = 


I  N  .  d  ■  .  , 
i  =  l  ^  ^ 


i  =  1,4 


( C . 14 ) 


and  performing  an  integration  by  parts  on  the  second  order 
derivatives  yields, 


//  NR,  dA  =  (l-p)N(k  )  (T  )  dl  +  //  ( 1-p)  (k  )  N  (T  )  dA 
A  "  xc  Ze  -ecn  A  e~rcr 

e 


~  //  ii^-N(ke)dA  +  //  (1-p)  (ke)Nz(Tc)zdA 


-j>  (1-p)  (k  )N(T  )  dA  +  //  hZN  (T  -T  )  dA 

ie  e  '  c  n  A  -  c  a 

e 


//  R  ZNdA  +  //  p  p  N  T  dA 
a  9  a  c  ~  c 


(C.15) 
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where  the  coefficients  of  the  response  variables  are  them¬ 
selves  functions  of  the  response  variables,  and  thus,  the 
equations  are  nonlinear.  In  accordance  with  the  Galerkin 
method,  the  final  system  of  ordinary  differential  equations 


was  obtained  by  setting  each  residual, 


each  basis  function, 


that  is , 


,  orthogonal  to 


/  /  N  Ri 

A 

e 


dA  = 


(C. 12 ) 


The  4 *NSNP  ordinary  differential  equations  given  by  Equations 
C.12  retain  the  character  of  the  original  set  of  partial 
differential  equations,  i.e.,  self-adjoint  operators  yield 
symmetric  matrices  and  non  self-adjoint  operators  yield 
nonsymmetric  matrices.  Thus,  linear  field  operators  trans¬ 
form  to  matrix  operators  and  nonlinear,  coupled  algebraic 
operators.  Incorporation  of  the  boundary  conditions  resulted 
in  4 *NSNP  nonlinear  coupled  ordinary  differential  equations, 


F  (ip,ip ,  t)  =  A  +  B*  -  F  + 


;c.i3) 


subject  to  initial  conditions,  where  A  is  a  ( 4 *NSNP ) * (4 *NSNP) 

matrix,  B  is  the  matrix  associated  with  the  linear  field 

operators  in  Expression  C.5,  F  is  an  excitation  vector, 

★ 

is  a  3  x  9  matrix  arising  from  a  bilinear  operator  treat¬ 
ment  of  the  reaction  terms  in  Expressions  C.8  and  C.ll. 
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NSNP 

.  (n  , e  ,  t )  =  l  N  .  (n  /  e )  9 ,  •  ( t) 

J  i=l  3  13 


(C.l) 


T  =  1  0  .  ( n  ,  £  ,  t ) 


NSNP 

l  N.  (n,E)  0O  .  (t) 


p  =  v  3  j  ( n  ,  e  » t  ] 


NSNP 

[  n  (n,E)0,.(t) 
j=l  3 


14-  ( n  r  E  ,  t ) 


NSNP 

[  N.  (n,e)  04.(t) 
j=l  D  43 


(C.  4) 


where  N^  for  j  =  1,...,NSNP  is  a  set  of  specified  linear  basis 
functions  with  local  support,  and  the  sets  6^,  e2j'  93j' 
and  1^.;  j  =  1 , .  .  . , NSNP ,  are  the  solution  coefficients  to  be 
determined.  The  N^  were  selected  to  satisfy  the  condition 
N j ( NP i )  =  5^^  where  the  Kronecker  delta,  6  ^  ^ ,  is  defined  by 
=  1  for  i  =  j,  and  5^  =  0  for  i  ^  j.  As  a  result,  6^, 
~2j'  ^3j  and  94j  are  t^ie  values  i>2'  ^3  and  ^4  at  the  nodal 

points  (i.e.,  ^  i  j  ( n  ,  e  ,  t )  =  Q^U)). 

Area  interpolation  functions  (shown  in  Figure  C.l) 
were  used  as  the  linear  basis  functions  which  provide  the 
necessary  function  continuity.  As  a  measure  of  error,  a 
residual,  R^,  is  defined  for  each  field  equation  by, 


R.  =  'p  +  A  .  (<p)  -  F  +  C  [T_4>]  * 

-L  ~  1  ~  ~  ^ 


( C  .  5 ) 
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from  the  system  matrices,  A(t),  B(t)  and  C(t).  In  order 

to  include  effects  of  T  and  f  arising  from  reaction  rate 

c 

terms  in  the  Jacobian  vector,  and  thereby  improve  the  effi¬ 
ciency  of  the  integration  routine,  the  combustion  terms  are 
incorporated  into  the  residual  equations  (in  DIFMOD) . 
Modifications  to  the  Jacobian  matrix  are  accomplished  in  the 
JACMOD  and  NUITSL  routines.  Reaction  rate  terms  of  Expres¬ 
sions  C.40  and  C.41  contribute  nine  terms  to  the  respective 
residual  equations  and  twenty-seven  terms  to  the  Jacobian 
vector.  Reaction  terms  are  generally  expressed  as, 

3  3 

RT  =  l  l  C..  ,  k  =  1,9  (C.53) 

i=!  j=l  1K  11  43 


3RT 

3'iu 


3 

Cik*  '^4 j 

3  =  1 


(C. 55) 


where  the  k*  are  compatible  with  the  column  locations  in  the 
C  matrix  of  the  products,  and 
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( C . 56 ) 


where  the  k*  are  compatible  with  the  column  locations  in  the 
C  matrix  of  the  products.  Terms  from  Equation  C.55 

may  be  incorporated  into  the  PW  Jacobian  vector  (containing 
contributions  to  the  Jacobian  from  the  linear  spatial 
operators)  as  contributions  from  Equation  C.S  with  combustion. 
Similarly,  terms  from  Equation  C.56  may  be  incorporated  into 
PW  for  terms  arising  from  the  0^  residual  equations  with 
combustion.  The  remaining  contributions  to  the  Jacobian 
vector,  Equation  C.56  for  the  carbon  equations  and  terms 
arising  from  Equation  C.55  for  the  0 ^  equations  are  stored 
in  the  PWMOD  (NDOF/2,9)  matrix.  The  column  number  or  variable 
of  differentiation  array,  INMOD  (NDOF/2,9)  and  PWMOD  array 
are  communicated  to  the  NUITSL  routine  so  additional  Jacobian 
terms  not  assimilated  into  the  PW  array  (by  virtue  of  the 
storage  scheme  selected)  may  be  taken  into  account  during 
the  convergence  sequence  for  the  new  iterates.  The  iterating 
scheme  is  of  Newton-Raphson  type.  Thus,  the  final  synthesis 
of  the  Jacobian  including  effects  of  all  terms  arising  from 
combustion  is  consummated. 

5 .  Derivation  of  the  FEM  Operators 

In  the  section  on  the  finite  element  formulation, 
the  following  six  differential  operators  were  identified. 
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if  N  (  ■+■)  dA 

A  ^ 


(c. 57; 


if  N (v)  dA 


(C. 58: 


//  N  ('{,■  )  dA 

A  ~r  r 


[C .  59 ) 


//  N  (’Jj)  dA 
A  ~  2  2 


;c. 60) 


/  N  ( -^ )  dA  and  J  j  N  ^  dA 


(C .  6 1 ) 


£  kN(^)  dH 

Ze  ~ 


(C . 62 ) 


where  N ^  are  the  global  basis  functions.  These  operators 
are  constructed  on  the  element  level  by  introducing  the 
corresponding  element  basis  functions,  The  global  and 

element  basis  functions  are  related  by, 


T  T 

■^i  =  N  0  +  H 


(C . 63 ) 


The  derivation  of  the  local  elemental  matrices  (using 
the  local  coordinate  system  depicted  in  Figure  C.3)  according 


to  the  Galerkin  method  for  the  global  operations  proceeds 
as  follows: 


For  Operator  C.29(also  C.57), 


Global 


Local 


//  N  (•.!,)  dA 
A  “ 


11  5  n  dA 

A  '  ~r~ 


(C. 64: 


Noting  that, 


c  .  9  . 


;  j  =  l ,  3 


( C  .  6  5  ) 


where  the  repeated  index  implies  Einsteinian  notation,  the 
elemental  matrix  becomes, 


//  5  J-edA 

A  e 


(C. 66! 


Expanding , 


b  .  T 


( 2A~)  ~9dA 
e 


bl  b2  b3 

k  bl  b2  b3 

bl  b2  b3 


e  (c.67; 


For  Operator  C.30  (also  C.58), 


Global 


Local 


//  N  (ij> )  dA 


//  5  5_0dA 

Ae  ~  Z~ 


(C. 68: 
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Substituting  the  local  shape  functions  gives, 


C .  T 

//  §<3aL>  5 dA 


A  e 

e  - 


;c. 69) 


the  elemental  matrix  becomes, 


C .  T 


e 


C  C  C 

C1  2  ^3 

6  C1  C2  C3 

C1  C2  C3 


(C. 70) 


For  Operator  C. 31 (also  C.59) 


Global 


Local 


//  N  (*)  dA 
A  ~r  r  e 


//  dA 

Ae 


(C . 71 ) 


Substituting  the  local  shape  functions  gives, 


b  .  b. 

^  ^  2A1"  2A1  ?  dA 

e  e 

A  . 

e 


4Aq 
e  e 


blb2  blb3 


772  V  blb2  b2 


b2b3 


blb3  b3bl 


( C .  7  2  ) 


the  elemental  matrix  becomes. 


If  N  ( !|> )  dA  ~ 
A  ' r  r 


blb2  blb3 


blb2 


b2b3 


blb3  b2b3 


(C  .  7  3  ) 
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For  Operator  C.32  (also  C.60), 


Global 


Local 


//  N  ( --j )  dA 

A  - 2  2 

e 


il  ?  dA 


( C  .  7  4  ) 


Substituting  the  local  shape  functions  gives. 


Ae 


cT 


2A 


S  dA  = 


II 


4A  A 
e  e 


C1  C1C2  C1C3 
C1C2  C2  C2C3 
C1C3  C2C3  C3 


edA  (C .  75 ) 


the  elemental  matrix  becomes. 


[I  h  " 

e 


1 

4A 


C1C2 

C1C3 


C1C2 


C2C3 


C1C3 

C2C3 


(C. 76) 


For  Operator  C.33  (also  C.61) 


Global 


Local 


II  N  ij;  dA 
A  ~ 


//  £  5  9  dA 

A  ~  " 


(C . 77 ) 


Substituting  the  local  shape  functions  gives, 
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// 


0  dA  = 


// 

A 


H  ’  2 


>2^3 


n’3 

f  ?• 
^2-3 


0  dA 


(C. 78) 


the  elemental  matrix  becomes, 


//  ST  0  dA  =  ^ 


2  11 
12  1 
112 


(C. 79) 


The  last  operator  C.34  (also  C.62)  is  incorporated  into  the 
excitation  vector  as  described  in  the  FEM  formulation  and 
has  been  addressed  in  Subsection  2,  Implementation  of 
Boundary  Conditions. 
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